Final  Performance  Report 


Contract  #  FA9550-08-C-0006: 

Hearing  Protection  for  High-Noise  Environments 
Period  of  Performance:  Oct  01,  2007  —  Nov  30,  2009 


Attachment  4 

Parallelization  of  the  Acoustic  Integral-Equation  Solver 
and  Example  Applications 


Prepared  by: 

MONOPOLE  RESEARCH 
739  Calle  Sequoia,  Thousand  Oaks,  CA  91360 
tel:  (805)  375-0318  fax:  (805)  499-9878 


Approved  for  public  release;  distribution  unlimited 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2oq<)  2.  REPORT  TYPE 

4.  TITLE  AND  SUBTITLE 

Parallelization  of  the  Acoustic  Integral-Equation  Solver  and  Example 
Applications 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

MONOPOLE  RESEARCH, 739  Calle  Sequoia, Thousand  Oaks, CA, 91360 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

41 

3.  DATES  COVERED 

01-10-2007  to  30-11-2009 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Contents 


1  Introduction  1 

2  The  FFT-based  matrix  compression  scheme  4 

2.1  Problem  discretization  and  iterative  solution  .  4 

2.2  Compressed  stiffness  matrix  representation  .  5 

2.3  Fast  matrix- vector  multiplication .  6 

3  Parallel  algorithm  implementation  8 

3.1  Data  partition  .  8 

3.1.1  Distribution  of  the  Cartesian  vectors .  8 

3.1.2  Distribution  of  the  geometry  data  .  10 

3.1.3  Distribution  of  the  MoM  unknowns .  11 

3.1.4  Distribution  of  the  near- field  matrix .  11 

3.1.5  Distribution  of  the  far-field  matrix  data .  11 

3.2  Distribution  of  the  computational  work  .  12 

3.2.1  Geometry  processing .  13 

3.2.2  Matrix  construction .  13 

3.2.3  Matrix-vector  multiplication .  15 

4  Computational  complexity  estimates  17 

5  Storage  estimates  17 

6  Examples  of  large-scale  computations  22 

6.1  The  geometrical  models  and  material  properties .  22 

6.2  Analysis  of  the  accuracy  of  the  far-field  expansion  in  the  problem  of  the 

human  head  helmet  with  cork  filling .  24 

6.3  Pressure  distributions  in  a  human  head  model:  dependence  on  the  discretization  28 

6.4  Pressure  distributions  in  a  human  head  model  in  the  presence  of  an  unattached 

helmet .  28 

6.5  Pressure  distributions  in  a  human  head  model  in  the  presence  and  absence 

of  a  helmet  with  a  cork  lining .  29 

6.6  Scaling  in  parallel  computations .  31 

A  Handling  of  zero  padding  and  communication  in  parallel  implementation 
of  FFTs  34 


References 


38 


List  of  Figures 


1  Partition  of  the  Cartesian  grid  and  the  object  geometry  into  disjoint  z-slices  10 

2  Model  of  the  external  head  surface  and  the  helmet .  22 

3  Averaged  (over  distance  bins)  relative  far-held  errors  in  the  surface  (S)  and 
volume  (V)  problems  for  the  range  parameter  RN/h  =  12,  plotted  as  func¬ 
tions  of  the  distance  between  basis  function  supports,  measured  in  units  of 
the  Cartesian  grid  spacing  h.  The  plots  represent  output  of  the  processor 

p  =  59  with  the  assigned  N ^  ~  60,  000  unknowns .  25 

4  Convergence  (the  relative  residual  norm)  of  iterative  solutions  in  the  surface 

(S)  and  volume  (V)  problems  for  the  indicated  values  of  the  near-held  range 
parameter  RN/h .  26 

5  Distributions  of  the  real  part  of  pressure  on  the  external  surface  of  the  model 
of  the  human  head,  the  helmet,  and  the  cork  filling,  with  N  ~  4,  700, 000 
unknowns.  The  distributions,  plotted  in  linear  scale,  were  obtained  for  the 
same  values  of  the  near-held  range  as  in  Fig.  4,  i.e. ,  (a)  R^/h  =  6,  (b)  R^/h  = 

8,  (c)  RN/h  =  10,  (d)  RN/h  =  12 .  26 

6  Distributions  of  the  absolute  value  of  pressure  in  the  coronal  section  of 
the  model  of  the  human  head,  the  helmet,  and  the  cork  filling,  with  N  ~ 

4,  700,  000  unknowns.  The  distributions,  plotted  in  linear  scale,  were  com¬ 
puted  for  the  same  values  of  the  near-held  range  as  in  the  previous  Figures, 

i.e.,  (a)  RN/h  =  6,  (b)  RN/h  =  8,  (c)  RN/h  =  10,  (d)  RN/h  =  12 .  27 

7  Distributions  of  the  real  part  of  pressure  in  the  coronal  section  of  the  model 
of  the  human  head,  the  helmet,  and  the  cork  filling,  with  N  ~  4,  700, 000 
unknowns.  The  distributions,  plotted  in  logarithmic  scale,  were  obtained 
for  for  the  same  values  of  the  near-held  range  before,  i.e.,  (a)  R^/h  =  6, 

(b)  RN/h  =  8,  (c)  RN/h  =  10,  (d)  RN/h  =  12 .  27 

8  Distributions  of  the  real  part  of  the  pressure  on  several  section  in  the  axial 

plane  of  the  human  head  model .  28 

9  Distribution  of  the  real  part  of  the  pressure  on  several  sections  in  the  axial 

plane  of  the  human  head  model  with  an  unattached  helmet .  29 

10  Pressure  distributions  in  the  coronal  plane  for  (i)  the  human  head  model  and 

(ii)  the  system  consisting  of  the  human  head  and  a  steel  helmet  models  .  .  31 

11  Matrix- vector  multiplication  time  as  a  function  of  the  number  of  Cartesian 

grid  nodes .  33 

12  FFT  in  the  (x,y)  directions  of  the  distributed  data  (a),  followed  by  transpo¬ 
sition  and  FFT  in  the  z  direction  (b) .  34 

13  Original  FFT  data,  their  rearranged  storage,  and  the  result  of  their  transpo¬ 
sition  by  means  of  the  MPI  all-to-all  communication  routine .  37 


ii 


1  Introduction 


We  have  recently  developed  a  fast  iterative  volumetric  integral-equation  solver  for  acoustic 
scattering  and  propagation  problems  involving  inhomogeneous  media.  The  solver  incor¬ 
porates  a  Fast  Fourier  Transform  (FFT)-based  matrix  compression  technique  (AIM)  and 
the  corresponding  fast  matrix- vector  multiplication  algorithm  [1,  2].  The  method  consists, 
essentially,  of  the  method-of-moments  (MoM)  computation  of  the  near-field  interactions 
involving  physical  sources  and  fields,  and  evaluation  of  the  far-held  by  means  of  the  FFTs 
operating  on  equivalent  source  and  held  distributions  dehned  on  nodes  of  a  regular  Cartesian 
grid.  The  FFT-based  matrix  compression  is  particularly  well  suited  to  volumetric  discretiza¬ 
tion  (with  a  fairly  uniform  discretization)  and  to  sub-wavelength  problems,  i.e. ,  spatial  res¬ 
olution  scale  much  smaller  than  the  wavelength.  Both  of  these  features  are  characteristic 
of  realistic  anatomical  models  and  propagation  problems  in  the  biological  applications  we 
are  considering. 

For  a  problem  with  N  unknowns  (in  this  implementation,  N  tetrahedra  in  the  volumetric 
mesh)  the  total  storage  is  O(N)  and  the  solution  complexity  as  0(N  log  N)  times  the 
number  of  iterations  (which  may  be  substantial,  but,  in  practice,  several  orders  of  magnitude 
smaller  than  the  number  of  unknowns).  Thus,  scaling  of  our  solution  method  is  similar  to 
that  of  differential-equation  (finite-element)  solvers.  We  have  also  carefully  optimized  the 
choice  of  the  matrix  compression  parameters,  in  order  to  minimize  the  total  cost  of  near- 
and  far- field  computations. 

Nevertheless,  in  application  to  large  problems  (millions  of  unknowns),  parallelization 
of  the  solver  is  essential.  We  describe  here  such  a  parallel  implementation  developed  for 
distributed-memory  systems  based  on  the  message  passing  interface  (MPI). 

The  most  essential  part  of  the  matrix  compression  and  fast  matrix-vector  multiplication 
procedures  is  the  FFT  algorithm.  Its  distributed-memory  parallelization  is  nontrivial,  since 
the  nature  of  the  algorithm  gives  rise  a  large  amount  of  inter-processor  communication:  even 
if  the  data  are  initially  partitioned  across  the  processors,  the  fact  that  Fourier  transforms 
have  to  be  evaluated  in  all  the  spatial  coordinate  directions  requires  a  global  rearrangement 
(“transposition”)  of  the  data.  Computation  of  the  near-field  contributions,  on  the  other 
hand,  is  to  large  extent  local,  and  requires  only  a  modest  amount  of  communication. 

For  these  reasons  we  build  the  parallelization  scheme  of  the  code  around  the  parallel 
FFT  algorithm.  We  take  here  advantage  of  the  availability  of  a  widely  used  FFT  package, 
FFTW,  which  allows  operations  on  FFT  data  distributed  spatially  (in  the  form  of  “slices”) 
across  the  processors.  Its  MPI-parallelized  implementation  is  available  in  both  the  current 
version  3  and  the  previous  version  2.  We  opted  for  the  more  recent  MPI  alpha  version 
3.2alpha3  (recently,  in  November  2008,  this  version  has  been  replaced  by  3.3alphal). 

Most  of  the  MPI-based  operations  are  transparent  to  the  package  user.  However,  a 
literal  application  of  the  FFTW  routines  to  carry  out  a  global  distributed  FFT,  would  have 
been,  in  our  case,  rather  inefficient.  The  main  reason  is  as  follows: 

(i)  The  far-held  computation  amounts  simply  to  a  convolution  of  the  given  distribution 
of  sources,  dehned  on  a  Cartesian  grid,  (we  describe  them  by  a  vector  X )  with  a  Green 
function  G.  This  convolution  can  be  evaluated  in  terms  of  FFTs.  However,  in  order  to  avoid 
aliasing  due  to  periodicity  of  the  FFT  (as  a  discrete  Fourier  transform),  the  “physical”  region 
on  which  X  is  dehned  has  to  be  enlarged  by  the  factor  of  two  in  all  directions,  and  padded 
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with  zeros. 

(ii)  In  a  parallel  implementation  the  padded  data  X  are  partitioned  into  slices  (say,  in 
the  direction  of  the  z  axis)  and  distributed  one  slice  per  processor.  Hence,  about  one-half 
of  the  processors  handles  only  the  padding  region. 

(iii)  Such  a  distribution  of  data  does  not  impair  in  any  way  execution  of  the  FFTs. 
However,  the  vector  multiplication  by  the  compressed  matrix  includes  also  near-field  oper¬ 
ations  which  involve  the  physical  sources  (say,  x)  and  field,  which  are  coupled  only  to  the 
“physical”  part  of  the  vector  X,  which  resides  only  on  one-half  of  the  processors.  Those 
near-field  computations  involve  matrix  blocks  which  have  substantial  sizes  (especially  in 
a  volumetric  problem).  One  possibility  is  to  store  those  blocks  only  on  that  half  of  pro¬ 
cessors  which  handle  the  physical  part  of  the  vector  A;  in  this  case  those  processors  are 
overloaded  with  storage,  while  the  memory  resources  of  the  remaining  processors  are  not 
fully  used.  Another  possibility  is  to  reorganize  the  near-field  matrix  blocks  and  distribute 
them  evenly  across  the  processors;  this  scheme,  however,  requires  an  additional  (and  com¬ 
putationally  expensive)  communication  between  practically  all  processors  before  and  after 
near-field  computation. 

In  view  of  the  above  difficulties,  we  devised  an  alternative  distribution  of  FFT  data, 
in  which  each  processor  stores  equal  parts  of  the  physical  part  of  the  Cartesian  vector  X 
and  the  padding  region.  With  this  “interleaved”  storage,  the  physical  sources  x  and  near¬ 
field  matrix  blocks  can  also  be  uniformly  distributed  across  the  processors,  and  no  extra 
communication  costs  are  incurred.  Additional  advantages  of  the  “interleaved”  data  layout  is 
a  reduction,  by  a  factor  of  about  two,  in  the  amount  of  communication  in  the  parallel  FFT 
(this  aspect  is  discussed  in  Appendix  A)  and  the  fact  that  it  allows  us  to  directly  utilize  the 
collective  MPI  communication  routine  MPI_Alltoallv,  which  is  often  highly  optimized  for 
specific  interconnect  networks. 

We  also  have  to  note  here  two  potential  drawbacks  of  the  interleaved  storage: 

(a)  Each  processor  has  to  store  not  only  the  “physical”  part  of  its  Cartesian  grid,  but  also 
the  padding  area.  In  the  alternative  scheme  (with  separate  physical  grid  and  padding 
storage)  the  full  (padded)  grid  can  be  distributed  into  twice  as  many  processors. 

(b)  The  FFTs  cannot  be  computed  entirely  in-place,  and  hence  an  additional  communication- 
related  buffer  is  needed. 

We  address  these  problems  in  Section  5  and  show  that,  in  practice,  they  do  not  lead  to 
difficulties. 

Scaling  of  a  parallel  implementation  of  the  solver  for  large  problems  requires  distributing 
of  practically  all  sizable  data  across  the  processors,  with  a  minimal,  if  any,  replication; 
at  the  same  time,  processors  must  have  access  to  some  non-local  geometry  data  in  order 
to  evaluate  near-field  couplings  between  sources  and  field  variables  assigned  to  different 
processors.  We  solve  the  problem  by  (i)  partitioning  the  geometry  into  non-overlapping 
slices  corresponding  to  the  Cartesian  grid  slices  assigned  to  the  processors,  but  (ii)  in  the 
matrix  construction  stage  temporarily  combining  (“welding”),  on  each  processor,  several 
slices  into  a  single  “stack”  geometry  of  thickness  equal  at  least  to  the  near-field  range.  This 
“welding”  procedure,  which  can  be  realized  in  a  straightforward  way  in  the  “interleaved” 
storage  approach,  allows  us  to  treat  the  entire  stack  of  slices  as  a  conventional  bona  fide 
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complete  geometry,  without  having  to  introduce  any  additional  connectivity  data  relating 
adjacent  or  identical  geometry  elements  in  different  slices  (which  would  have  been  necessary 
if  slices  were  treated  as  separate  geometries). 

We  describe  here  the  details  of  parallelization  in  the  present  code,  which  assumes  basis 
functions  associated  with  tetrahedra.  However,  these  parallelization  techniques  are  applica¬ 
ble,  with  only  minor  modifications,  to  various  other  types  of  geometry  discretization,  e.g., 
to  basis  functions  associated  with  nodes  of  the  mesh. 
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2  The  FFT-based  matrix  compression  scheme 

We  give  here  only  a  brief  summary  of  the  FFT-based  compression  method  used  in  our 
solver.  More  details  can  be  found  in  Refs.  [2]  and  [3]. 

2.1  Problem  discretization  and  iterative  solution 

We  solve,  iteratively  (by  using  one  of  the  minimum-residual  algorithms,  GMRES  or  GCR 
[4,  5]),  the  linear  system 

Ax  =  b  (2.1) 

obtained  by  Galerkin  discretization  of  the  acoustic  Lippmann-Schwinger  volumetric  integral 
equation  for  the  pressure  p , 

Kp  =  p inc,  (2.2) 

describing  scattering  of  an  incident  pressure  wave  pmc  on  a  certain  bounded  region  fl  filled 
with  an  inhomogeneous  material  characterized  by  a  density  p(r)  and  the  compressibility 
«(r).  According  to  the  nature  of  the  L-S  equations,  the  Green  function  appearing  in  the 
integral  operators  is  always  that  of  the  background  (surrounding)  space  -  here  air.  Corre¬ 
spondingly,  the  equations  involve  ratios  of  the  scatterer  density  and  compressibility  to  those 
of  the  background  medium. 

In  fact,  in  the  case  of  large  density  contrasts,  we  use  a  modified  integral-equation  formu¬ 
lation  described  in  Ref.  [3].  It  involves  three  different  operators,  associated  with  a  surface 
problem  (defined  on  high-contrast  surfaces),  with  a  coupling  between  that  surface  and  the 
volume  region,  and  with  a  volume  problem.  Explicit  expressions  for  these  operators  are 
also  given  in  Ref.  [2]. 

The  vector  x  consists  of  coefficient  of  expansion  of  the  pressure  p  into  a  set  of  local  basis 
functions  cj)al  and  b  is  a  vector  obtained  by  projecting  the  incident  field  pmc  on  the  same  set 
of  basis  functions.  The  elements  of  the  stiffness  matrix  A  are  given  by  the  MoM  expression 

Aap  =  (0a  ,  K  0p)  =  J  d3?T  d3r2  ^>a(r1)  K(rv  r2)  ^(r2)  ,  (2.3) 

where  K  is  one  of  the  (non-symmetric)  kernels  arising  in  our  system  of  the  integral  equations. 
Similarly,  the  elements  of  the  r.h.s.  vector  b  are 

K  =  (K,  P'"c )  =  / d3r  <pa(r) p”‘(r)  .  (2.4) 

The  present  implementation  of  the  code  assumes  constant  basis  functions  supported 
on  tetrahedra,  hence  the  indices  a  and  f3  are  tetrahedron  numbers.1  We  assume  that  the 
material  properties  of  the  medium  (its  density  p  and  compressibility  n)  are  constant  on  the 
individual  tetrahedra. 

lrThe  general  structure  of  the  code  and  the  compression  algorithm  remains  unchanged  for  other  types  of 
basis  functions,  e.g.,  piecewise-linear  functions  associated  with  vertices. 
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2.2  Compressed  stiffness  matrix  representation 

Since,  in  an  iterative  solver,  most  operations  amount  to  matrix-vector  multiplication,  a  fast 
solution  can  be  achieved  through  matrix  compression  and  a  fast  matrix- vector  multiplication 
algorithm. 

The  FFT-compressed  matrix  A  is  stored  as 

Aal ,  A$T  +  .43?  .  (2.5) 

The  near-field  part  of  the  matrix,  coupling  tetrahedra  up  to  the  relative  distance  of  some 
near-field  range  i?N,  is  computed  by  using  the  conventional  MoM,  i.e.,  Eq.(2.3)  and, 
according  to  the  usual  AIM  prescription  [1,  2],  subtracting  the  far- field  contributions  (as 
defined  below).  At  larger  distances  the  couplings  are  described  by  the  far- field  part  of  the 
matrix,  having  the  general  structure 

A“  =  E  r™  G(u  -  v)  y£> ,  (2,6) 

CSS 

where  u  and  v  are  nodes  of  a  global  Cartesian  grid  C  enclosing  the  object,  and  B{a)  is 
a  cubic  set  of  ( M  +  l)3  grid  nodes  at  which  equivalent  sources  associated  with  the  basis 
function  a  are  located;  we  call  this  set  the  expansion  box  of  the  basis  function  a  (and 
similarly  for  B((3)).  The  parameter  M  is  the  multipole  expansion  order  of  the  far  field; 
the  typical  value  M  =  2  corresponds  to  octupole  expansion.  For  each  tetrahedron  a  its 
expansion  box  B(a)  is  defined  as  that  set  of  (M  +  1)  x  (M  +  1)  x  (Af  +  1)  grid  nodes,  whose 
center  is  closest  to  the  centroid  of  the  tetrahedron. 

The  expansion  coefficients  and  vftj  relate  the  original  MoM  basis  functions  and 
sets  of  the  equivalent  sources,  and  also  depend  on  the  structure  of  the  integral  operator  in 
question.  These  coefficients  are  determined  by  requiring  that  the  original  basis  functions 
and  the  corresponding  sets  of  equivalent  sources  have  the  same  multipole  moments  up  to 
some  expansion  order  M.  Finally,  G(u  —  v)  are  the  tabulated  values  of  the  (scalar)  Green 
function  of  the  Helmholtz  equation  with  the  incident  field  wave  number  k  (describing  wave 
propagation  in  the  background  medium). 

Our  formulation  for  high-contrast  problems  involves  matrix  elements  of  two  types, 
“monopole”  (M)  and  “dipole”  (D),  for  which  the  representation  (2.6)  has  the  form 

43/"  =  k2  E  C  G(u  -  V)  <’> 

U  v  '  'ft  ^  ' 

A/"  =  -  E  '  G<u  -  v)  Vlv  • 

11, V 

(Eqs.  (21)  and  (24)  of  Ref.  [2]).  It  involves  three  different  types  of  expansion  coefficients, 
with  one  of  them,  V^,  depending  on  the  discontinuities  A  of  the  density  across  facets  of  the 
tetrahedron  (5\  similarly,  and  Kp  are  the  density  and  compressibility  of  the  tetrahedron 
f3,  and  p0  and  n0  the  corresponding  parameter  values  in  the  background  space  (air).  The 
vector- valued  expansion  coefficients,  and  V^,  represent  the  gradient  of  the  pressure. 


(2.7a) 

(2.7b) 
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Having  specified  the  far-field  representations  (2.6)  or  (2.7)  of  the  matrix  elements,  we 
define  the  near-field  component  of  the  matrix  as  a  difference  between  the  actual  (MoM) 
matrix  elements  and  their  far-field  approximations, 


A 


Near 

aj3 


4*/3  ^a/3  f°r  lc«  c/3 1  — 

0  otherwise  , 


(2.8) 


where  cQ  and  Cg  are  centroids  of  the  supports  of  the  basis  functions  <fia  and  4>p,  and  i?N  is 
a  pre-defined  near-field  range.  Subtraction  of  the  far-field  matrix  element  A^g  (mentioned 
following  Eq.(2.5))  is  necessary  in  view  of  the  fact  that,  in  matrix-vector  multiplication,  the 
far-field  contributions  are  being  added  irrespectively  of  the  separation  between  the  basis 
functions  a  and  (3  (this  is  required  by  the  Toeplitz  property  of  the  compressed  matrix)  and, 
clearly,  AFgr  is  not  a  good  approximation  to  Aap,  hence  the  subtraction  compensates  for 
the  resulting  error. 

According  to  Eq.(2.8),  for  distances  exceeding  /?,N,  we  neglect  the  difference  between  the 
matrix  element  Aan  and  its  far-field  approximation  AFgr.  Therefore,  the  near- field  range 
has  to  be  large  enough  to  ensure  an  acceptably  small  relative  error  of  the  approximation, 


^ ad 


I  A  _  4Far| 
\AaP  Aa/3\ 
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a/3 1 


(2.9) 


for  \ca  —  cJ  >  i?N.  Usually,  an  error  of  1  %  to  2 %  can  be  considered  adequate. 

In  the  following  we  denote  the  average  tetrahedron  size  (edge  length)  by  a,  and  the 
spacing  of  the  Cartesian  grid  by  h.  Our  experience  with  optimizing  compression  parameters 
(the  main  results  of  which  were  reported  in  Ref.  [2] )  has  shown  that  preferred  values  of  the 
grid  spacing  is  relatively  small,  h<a/ 2,  the  optimal  expansion  order  is  M  =  2,  and  the 
near-field  range  can  be  chosen  as  ~  3  h  for  the  monopole  matrix  elements  and  i?N  ~  6  h 
for  the  dipole  elements  (Eqs.  (2.7)).  With  these  parameters  the  relative  error  in  the  far-field 
matrix  elements  should  be  on  the  level  of  about  1  %. 

Eqs.  (2.5)  and  (2.6)  can  be  interpreted  as  providing  a  compressed  matrix  representation 
as  a  sum  of  the  near- field  matrix  (which,  by  construction,  is  sparse)  and  a  far-field  matrix 
given  as  a  product  of  two  sparse  matrices  V  and  a  Toeplitz  matrix  G.  For  N  unknowns  and 
Nc  ~  N,  and  for  a  fixed  near-field  range  i?N,  all  these  matrices  require  storage  proportional 
to  N ;  hence  the  entire  storage  is  also  O(N). 


2.3  Fast  matrix-vector  multiplication 

With  the  compressed  matrix  representation  of  Eqs.  (2.5)  and  (2.6),  matrix-vector  multipli¬ 
cation  y  =  Ax  amounts  to  computing 

y  =  y Near  +  y Far  =  ANear x  +  VGVTx  .  (2.10) 

While  the  near-field  evaluation  involves  simply  a  vector  multiplication  by  a  sparse  matrix, 
the  far-field  contribution  is  obtained  by 
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1.  Accumulating  contributions  of  basis  functions,  weighted  by  the  components  Xg  of  the 
vector  x ,  to  the  equivalent  sources  X  on  the  Cartesian  grid  nodes  v  e  C, 

Xv+=  Vjg  Xg  for  v  G  B{0)  .  (2.11) 

This  operation  is  carried  out  in  a  loop  over  (3,  and  in  each  step  contributions  to  all 
relevant  grid  nodes  v  are  computed  and  added.  As  indicated  in  Eq.(2.7),  the  vector 
X  and  the  coefficients  V  have,  in  our  case,  four  components,  describing  the  pressure 
and  its  gradient. 

2.  Taking  its  Fourier  transform 

E  =  E  -^qv  •  (2.12) 

vec 

3.  Multiplying  X  by  the  Fourier  transform  of  the  Green  function, 

yq  =  GqAq,  (2.13) 

for  all  points  q  of  the  dual  grid  C. 

4.  Taking  the  inverse  Fourier  transform 

=  E  -^uq1  Vq  ■  (2.14) 

qeC 

5.  Finally,  converting  back  the  Cartesian-grid  field  distribution  Y  to  its  MoM  represen¬ 
tation  y. 

y*™=  E  ^au^u-  (2-15) 

u  £B(a) 

Physically,  X  in  the  above  expressions  represents  the  equivalent  sources  on  the  Cartesian 
grid  C,  and  Y  is  the  field  generated  on  the  same  Cartesian  grid  by  the  sources  X.  We  recall 
that  the  vectors  X  and  Y  have  four  components,  while  G  is  a  scalar. 

In  the  following  we  refer  to  the  sets  of  equivalent  sources,  fields,  or  the  Green  function 
values,  and  to  their  discrete  Fourier  transforms,  i.e. ,  to  the  vectors  X,  Y,  G,  X,  Y,  or  G, 
as  Cartesian  vectors  (all  of  them  are  denoted  by  capital  letters).  We  note  that  in  our 
implementation  FFT  operations  are  executed  in-place,2  hence  the  vectors  X,  Y,  G,  X,  and 
Y  share  the  same  storage  area. 

2However,  as  we  mentioned  in  the  Introduction  and  discuss  later,  communication  in  the  parallel  FFT 
requires  an  additional  storage  buffer.  Its  size  is  one-half  of  the  vector  X,  since  zero-padding  in  one  direction 
is  not  needed. 
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3  Parallel  algorithm  implementation 


A  parallel  implementation  of  the  solver  for  large  problems  requires  distributing  of  practically 
all  data  across  the  processors,  with  a  minimal,  if  any,  replication.  Similarly,  the  compu¬ 
tational  work  should  be  distributed,  as  evenly  as  possible,  across  all  the  processors.  Our 
implementation  achieves  this  goal  and  ensures  scaling,  in  the  sense  that  no  processor  stores 
an  amount  of  data,  or  executes  an  amount  of  work,  growing  with  the  total  problem  size.  In 
other  words,  the  storage  and  computational  work  per  processor  may  remain  bounded  if  the 
number  of  processors  grows,  in  an  appropriate  way,  with  the  problem  size. 

We  discuss  in  the  following,  separately,  the  aspects  of  parallelization  related  to  data 
partition  and  distribution  of  the  computational  work. 

3.1  Data  partition 

The  data  to  be  distributed  across  the  processors  include: 

1.  the  object  geometry  Q, 

2.  the  MoM  unknowns  x  and  the  right-hand-side  vector  b, 

3.  the  equivalent  Cartesian  representations  of  the  unknown  vector  and  the  Green  func¬ 
tion,  i.e,  the  Cartesian  vectors  X,  Y,  X,  and  Y  (sharing  common  storage),  defined  on 
points  of  the  Cartesian  grid  C  covering  the  object, 

4.  the  elements  of  the  compressed  representation  (2.6)  of  the  stiffness  matrix  A: 

(a)  the  near-field  part  ANcar  of  the  matrix, 

(b)  the  expansion  coefficients  V  mapping  MoM  into  Cartesian  vector  representations, 

(c)  the  tabulated  Green  function  G,  or  rather  its  Fourier  transform  G  appearing  in 
Eq.(2.13). 

In  the  following  we  assume,  for  definiteness,  that  the  geometry  and  the  Cartesian  grid 
are  partitioned  into  “z-slices”  specified  by  the  ^-coordinate.  Below  we  describe  how  the 
data  elements,  listed  above,  are  distributed  across  a  selected  number,  P,  of  processors. 

3.1.1  Distribution  of  the  Cartesian  vectors 

As  we  mentioned  before,  the  matrix  compression  algorithm  and  its  parallel  version  are 
designed  around  the  operations  on  the  Cartesian  grid  data. 

The  global  Cartesian  grid.  Suppose  the  bounding  box  of  the  geometry  Q  has  sizes 
Bx  x  By  x  Bz ,  the  Cartesian  grid  spacing  (assumed  to  be  identical  in  all  three  directions) 
is  h,  and  the  far-held  expansion  order  is  M.  The  global  Cartesian  grid  C  is  then  defined 
as  the  set  of  Kx  x  Ky  x  Kz  grid  nodes  such  that  it  extends  outside  the  bounding  box  by 
at  least  (M  +  l)/2  grid  spacings  in  each  direction;  this  condition  is  necessary  in  order  to 
ensure  that  the  expansion  boxes  B(a)  of  all  the  tetrahedra  a  are  contained  in  the  grid  Q. 
The  physical  sizes  of  the  Cartesian  grid  are  then  (Kx  —  1)  h  x  (Kx  —  1)  h  x  (Kx  —  1)  h.  As 


we  discuss  in  Sections  4  and  5,  the  optimal  (from  the  point  of  view  of  the  Cartesian  data 
manipulations)  orientation  of  the  object  is  such  that 

I<x  <  Ky  <  Kz  ;  (3.1) 

we  always  assume  this  orientation  in  the  following. 

We  recall  that,  before  executing  the  Fourier  transform  of  Eq.(2.12),  the  Cartesian-grid 
data  X  have  to  be  extended  to  twice  their  “physical”  size  by  padding  them  with  zeros,  in 
order  to  avoid  aliasing  (due  to  periodicity  of  discrete  Fourier  transforms)  in  their  subsequent 
convolution  (2.13)  with  the  Green  function.  After  the  backward  Fourier  transformation  to 
the  coordinate  space  (Eq.(2.14)),  only  the  physical  part  of  the  output  vector  Y  is  used. 
With  padding,  the  numbers  of  Cartesian  grid  nodes  are  N.t  =  2  K i,  i  =  x,y,  z. 


Data  layout  in  the  parallel  FFT  algorithm.  Before  discussing  Cartesian  data  par¬ 
titioning  we  recall  that  the  FFTW  parallel  implementation  of  the  FFT,  say  the  Fourier 
transform  (2.12),  involves 

1.  Distributing  the  vector  X  across  the  processors  according  to  the  ^-coordinates  of  the 
grid  points. 

2.  Executing,  in  parallel,  the  two-dimensional  FFT  of  X  in  the  (x,  y)  directions. 

3.  Transposing  the  y  and  z  indices  in  the  data,  i.e.,  redistributing  the  data  such  that 
they  partitioned  according  to  their  y  indices.  This  operation  is  crucial  from  the  point 
of  view  of  the  parallel  algorithm  implementation,  since  it  involves  global  all-to-all 
communication  between  the  processors. 


4.  Executing,  in  parallel,  the  one-dinrensional  FFT  of  X  in  the  z  direction. 


After  the  transformation  the  data  are  arranged  in  a  different  than  the  initial  order,  but,  in 
our  application,  there  is  no  need  to  rearrange  them  (the  operations  on  the  Fourier  trans¬ 
formed  Cartesian  vectors  may  be  executed  on  the  reshuffled  data  as  well).  In  the  inverse 
FFT  the  steps  1  to  4  above  are  carried  out  in  the  reversed  order. 

An  important  remark  is  in  order  here.  We  mentioned  before  that,  in  our  acoustic  prob¬ 
lem,  the  Cartesian  vector  X  has  four  components,  describing  the  pressure  and  three  com¬ 
ponents  of  its  gradient.  However,  in  our  implementation  we  do  not  store  all  the  components 
of  that  vector  simultaneously,  but  rather  reuse  a  single-component  vector.  In  other  words, 
we  evaluate  the  far-held  matrix-vector  product  of  Eq.(2.10)  as  a  sequentially  computed  sum 
of  four  terms, 


Far 

a 


-  E  EECl’Wu-vivjy^, 

j=x,y,z  (3  u,v 


(3.2) 


where  we  used  Eqs.  (2.7)  and  labeled  with  the  index  j  the  vector  components  of  vj^  and 

Such  an  organization  of  the  data  results  in  a  very  significant  reduction  of  storage  and, 
in  practice,  no  loss  in  the  computation  time;  we  discuss  this  point  further  in  Section  3.2.3, 
4,  and  5. 
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Partition  of  the  Cartesian  grid.  According  to  the  requirements  of  the  parallel  FFT 
algorithm,  the  Cartesian  grid  Q  is  distributed  across  the  P  processors  by  partitioning  the 
set  of  K z  unpadded  grid  nodes  into  P  ^-slices,  i.e. ,  P  subsets  of  Kz  points,  such  that 

+  K +  •  •  •  +  Kf-V  =  Kz  .  (3.3) 


The  partition  task  is  carried  out  by  the  FFTW  routine 

f  f  tw_mpi_local_size_many_transposed.  The  same  routine  supplies  a  preferred  distribu¬ 
tion  of  the  Cartesian  vector  X  into  y-slices  used  in  the  z-direction  Fourier  transform  (step 
4  above),  i.e.,  a  partition 

+  IV«  +  ■  ■  •  +  lYf-1)  =  Ny  ,  (3.4) 

(p) 

in  this  case  specified  for  the  padded  grid  sizes,  in  practice,  most  of  the  numbers  Kz  and 
NyP ^  are  equal,  but  some  (those  assigned  to  the  last  processors)  may  be  zero. 

In  the  following  we  denote  by  C ^  the  z-slices  of  the  physical  Cartesian  grid  C  associated 
with  the  Cartesian  data  slices. 

3.1.2  Distribution  of  the  geometry  data 

Since  the  Cartesian  vector  data  are  obtained  from  the  MoM  data,  and  those  are  related  to 
the  object  geometry  Q.  we  partition  the  geometry  in  a  way  closely  related  to  the  Cartesian 
data  slices.  Specifically,  we  define  the  geometry  slice  Q  assigned  to  the  processor  p  is  defined 
as  the  collection  of  all  tetrahedra  whose  centroids  c  satisfy  the  condition 

zp  <  cz  <  zp+1  ,  (3.5) 

where  z0,  ...  ,zp  are  the  ^-coordinates  of  planes  located  half-way  between  the  nodes  of  the 
adjacent  grid  slices  C^p\  The  resulting  partition  is  illustrated  in  Fig.  1. 


Figure  1:  Partition  (for  P  =  3  processors)  of  the  Cartesian  grid  C  and  the  object  geometry  Q 
into  disjoint  z-slices  C ^  and  .  In  order  to  match  the  usual  vector  and  matrix  notation, 
we  directed  the  z  axis  downwards. 
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3.1.3  Distribution  of  the  MoM  unknowns 


Geometry  partition  determines  directly  the  distribution  of  the  vector  of  the  MoM  unknowns 
x.  It  is  imply  split  into  P  blocks, 


x  = 


x(°) 

^(i) 


(3.6) 


x 


OP-i) 


where  the  block  x (p)  corresponds  to  the  tetrahedra  in  the  geometry  slice  Q^p\ 


3.1.4  Distribution  of  the  near-field  matrix 

The  near-field  matrix  TNoar  (the  item  4(a)  on  the  list  in  Section  3.1)  must  include  couplings 
between  adjacent  geometry  slices,  up  to  the  near-field  range  f?N.  Thus,  generally,  klNear  has 
the  block  structure 


4(0,0) 

4(0,1)  . . 

4(0,0) 

0 

0 

4(ho) 

4(M)  . . 

4(1,^) 

4(1, 0+1) 

0 

4(0,0) 

4M)  . . 

4  (o’,' °0 

4(<T,<T+1) 

0 

(3.7) 

0 

4(0+1, 0)  . . 

•  4(f7+1'cr) 

4(o+l, o+l)  .  . 

0 

0 

0 

0 

0 

.  4(P-1,P-1) 

In  order  to  simplify  the  notation,  we  assumed  here  that  each  geometry  slice  couples  to  a 
of  its  nearest  neighbors,  i.e. ,  each  row  and  column  in  the  matrix  consists  of  up  to  2a  +  1 
blocks.  Since  slices  may  have,  in  general,  different  thickness,  the  number  of  nonzero  blocks 
in  tows  and  columns  may  vary.  In  the  following  a  will  denote  the  maximum  number  of 
coupled  slices.  We  note  that  in  our  problem  the  matrix  ylNear  is  not  symmetric. 

In  distributing  the  matrix  blocks  we  simply  assign  all  the  blocks  of  the  matrix  in  its  p- th 
row  to  the  processor  number  p  (assignment  of  blocks  by  columns  would  merely  amount  to 
changing  the  order  of  some  operations  in  matrix- vector  multiplication).  In  other  words,  the 
p-th  processor  stores  the  set  of  blocks 

J$Pi° 1  _  |4(P,P-o)  .  .  .  !)  J\(p,p)  j{(P,P+ 1)  .  .  .  ^(p.P+o jj  _ 

3.1.5  Distribution  of  the  far-field  matrix  data 

The  MoM-to-Cartesian  mapping  coefficients  V.  Partition  of  the  expansion  coef¬ 
ficients  V  (the  item  4(b)  listed  in  Section  3.1)  deserves  special  care,  since  their  storage 
requires  a  significant  amount  of  memory  (we  give  estimates  in  Section  5). 
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As  in  the  case  of  the  near-field  matrix,  the  coefficients  V au  may  couple  elements  assigned 
to  different  slices;  i.e. ,  a  tetrahedron  a  located  in  the  geometry  slice  Ql'p'1  may  contribute  to  a 
Cartesian  grid  node  u  belonging  to  a  grid  slice  C^p+T\  r/0.  Therefore,  a  given  processor  p 
must  store  some  V  data  associated  with  either  tetrahedra  or  Cartesian  grid  nodes  assigned 
to  other  processors. 

For  example,  a  processor  p  could  store  sets  of  coefficients  Vau  for  all  tetrahedra  a  such 
that  at  least  one  of  the  nodes  u  belongs  to  the  expansion  box  of  the  tetrahedron,  u  e  13(a) 
is  located  in  the  grid  slice  This  condition  implies  that  the  processor  may  have  to  store 
coefficients  V  for  tetrahedra  whose  centroids  are  up  to  the  distance  Mh  outside  the  home 
slice  of  the  processor.  Such  a  solution  might  be  acceptable  for  “thick”  slices  (of  thickness 
large  compared  to  h).  However,  for  thin  slices  (and  a  grid  slice  may  contain  just  a  single 
plane  of  the  grid),  the  additional  (and  replicated)  storage  could  easily  exceed  the  storage 
for  the  slice  itself. 

In  order  to  eliminate  replication  in  the  storage,  and,  at  the  same  time,  minimize  the 
amount  of  communication,  we  introduced  a  “sparse”  data  structure  for  the  representation  of 
the  mapping  coefficients  V :  It  includes,  for  each  tetrahedron  a,  a  list  of  nodes  u  (numbered 
within  the  expansion  box  13(a) )  for  which  the  coefficients  are  stored.  Thus,  each  processors 
stores  the  coefficients  for  all  tetrahedra  a  whose  expansion  boxes  overlap  the  the  grid  slice 
(j(p);  however,  it  stores  coefficients  only  for  the  nodes  u  €  C^p\  In  this  way,  expansion 
coefficients  for  each  a  and  u  are  stored  only  once,  and  the  additional  lists  amount  to  only 
a  very  minor  overhead  (more  detailed  estimates  are  given  in  Section  5). 

The  Green  function.  The  last  item,  4(c),  on  the  list  in  Section  3.1  is  the  Green  function 
tabulated  at  the  nodes  of  the  Cartesian  grid.  Its  multi-processor  distribution  is  obtained  in 
a  straightforward  way:  (i)  each  processor  p  computes  and  stores  data  G  for  its  own  grid  slice 
C<*\  and  (ii)  a  global  FFT  transform  is  applies  to  these  distributed  data.  The  result  is  the 
Fourier-transformed  g  distributed  according  to  the  y-slices  of  the  grid,  which  automatically 
matches  the  layout  of  the  Fourier  transform  X,  as  described  in  Section  3.1.1. 

3.2  Distribution  of  the  computational  work 

Parallelization  of  the  operations  in  the  code  is  constructed  around  the  multi-processor  par¬ 
tition  of  the  data.  An  important  feature  of  the  implementation  is  that  each  processor 
handles  only  the  minimal  part  of  the  object  geometry,  necessary  in  order  to  construct  the 
data  assigned  to  the  processor.  In  the  matrix  construction  stage  we  also  make  use,  rather 
freely,  of  disk  storage  of  geometry  data  and  the  matrix  blocks,  in  order  to  avoid  storing  too 
much  data  at  the  same  time.  In  the  iterative  solution  stage,  however,  we  do  not  utilize  any 
out-of-core  storage. 

In  addition  to  few  parameters  specifying  matrix  compression,  the  incident  wave,  etc.,  the 
input  data  consists  of  a  single  geometry  file  containing  vertex  coordinates  of  the  tetrahedral 
mesh  and  definitions  of  tetrahedra  as  quadruplets  of  vertex  numbers,  and  a  file  listing 
material  parameters  for  all  the  tetrahedra.  With  this  input,  the  computational  procedure 
involves  steps  described  in  the  following. 
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3.2.1  Geometry  processing 

The  main  operations  performed  here  are  as  follows: 

1 .  One  processor  extracts  from  the  input  geometry  file  the  minimal  amount  of  the  global 
geometry  information,  required  in  order  to  determine  geometry  and  other  data  parti¬ 
tion. 

2.  Each  processor  p  extracts  from  the  input  file  its  own  geometry  slice  and  stores 
it  as  a  separate  file.  It  then  extracts  the  corresponding  material  parameter  data  and 
also  stores  it  as  a  file. 

3.  Each  processor  p  reads  geometry  and  material  data  for  up  to  2<r  +  1  geometry  slices, 
g(p-v)  g(p+<r) ,  where  a  is  the  number  of  nearest  coupled  slices  (cf.  Eq.(3.7)).  In  the 
following  we  refer  to  the  slice  Q ^  originally  assigned  to  the  processor  p  as  its  “home 
slice” . 

The  processor  then  combines  the  set  of  slices  into  a  “welded  stack”  W’[p,<Tl.  We  use 
the  expression  “welded”  to  indicate  that  common  elements  of  adjacent  slices  (ver¬ 
tices  and  facets)  have  been  eliminated,  and  the  resulting  geometry  can  be  treated 
as  a  conventional  bona  fide  complete  geometry.  The  processor  also  concatenates  the 
corresponding  material  data  into  a  single  block. 

The  above  operations  provide,  for  each  processor,  the  input  data  needed  in  the  matrix 
computation. 

3.2.2  Matrix  construction 

The  near-field  matrix.  Each  processor  computes  then  the  near-field  part  of  the  matrix 
(Eq.(2.5))  by  using  a  pair  of  geometries,3  and  An  important  advantage  of 

this  procedure  is  that  any  matrix-fill  algorithm  developed  for  a  serial  code  may  be  applied 
here;  in  particular,  there  is  no  need  to  introduce  any  additional  connectivity  data  relating 
adjacent  or  identical  geometry  elements  in  different  slices.  The  resulting  rectangular  matrix 
block  constitutes  the  set  of  blocks  in  the  p-th  row  of  the  matrix  of  Eq.(3.7),  denoted  by 
am  in  Eq.(3.8).  However,  the  columns  of  this  block  matrix  are  still  indexed  by  tetrahedra 
in  the  welded  stack  geometry  yvdp,<Tl. 

Far-field  subtractions.  In  the  next  step,  each  processor  p  modifies  its  matrix  block  A^P'^ 
by  subtracting  far-field  contributions  evaluated  according  to  Eq.(2.6).  More  precisely,  the 
subtraction  is  implemented  as 

AV  AV  -  E  E  v£>G(u-r)v£>.  (3.9) 

ueZ3(a)  v£B(p) 


In  more  detail: 

3Here  and  in  similar  operations  the  geometries  are  always  complemented  with  their  corresponding  material 
data. 
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1.  We  tabulate  the  set  of  all  values  of  G(u)  needed  in  Eq.(3.9).  This  procedure  saves 
a  considerable  amount  of  time  by  eliminating  repeated  computation  of  exponents; 
the  required  storage  for  G  is  minor,  since  padding  is  not  needed  and  one  can  take 
advantage  of  the  symmetries  of  the  Green  function. 

2.  The  full  set  of  the  coefficients  Vau  is  computed  and  stored,  for  all  tetrahedra  a  in  the 
geometry  slice  G^  and  for  all  nodes  in  their  expansion  boxes  13(a). 

3.  A  loop  over  nonempty  columns  (3  of  the  matrix  block  A is  executed.  For  each 
tetrahedron  (3 : 

(<2\ 

(a)  The  set  of  coefficients  is  evaluated  and  stored  (in  a  small  work  array)  for  all 
nodes  v  in  the  expansion  box  B(/3). 

(b)  A  loop  is  executed  over  all  nonzero  elements  a  in  the  column  (3  of  the  matrix 
block  A^'g\  For  each  element,  associated  with  a  tetrahedron  a  in  the  geometry 
slice  G^p\  the  far- field  value  of  the  matrix  element  is  evaluated  according  to 

Eq.(3.9),  in  terms  of  the  previously  computed  and  stored  coefficients  and  the 

(2) 

computed  on-the-ffy  and  temporarily  stored  coefficients  ;  the  Green  function 
values  G(u  —  v)  are  extracted  from  the  table  precomputed  in  step  1  above. 

After  the  far-held  subtraction  is  completed,  the  set  of  mapping  coefficients  V^1)  and  the 
tabulated  G  are  deleted. 

The  advantage  of  the  above  procedure  is  that  it  requires  only  storing  the  mapping 

coefficients  for  tetrahedra  in  the  home  slice,  and  thus  avoids  any  replication  of  the 

(2) 

data  on  different  processors;  since  the  coefficients  VaJ  are  computed  on-the-ffy,  they  incur 
no  storage  whatsoever.  We  note  here  that,  for  thin  slices  (which  may  often  be  the  case), 

the  near-held  range  may  be  signihcantly  larger  than  the  slice  thickness;  hence,  storing  the 
(2) 

coefficients  for  (3s  in  the  entire  near-held  range  would  require  an  excessive  amount  of 
memory. 

Far-field  data  in  the  compressed  matrix.  Now  the  processor  p  continues  to  utilize  the 
availability  of  its  welded  stack  geometry  yV’[p,f7^  to  (re)conrpute4  the  hnal  sets  of  the  mapping 
coefficients  and  .  This  time  the  coefficients  are  stored  in  the  sparse  representation, 
as  dehned  in  Section  3.1.5;  i.e.,  they  are  computed  only  for  nodes  u  in  the  home  grid  slice 
C^p\  and  for  all  necessary  tetrahedra  a,  either  in  the  home  geometry  slice  or  in  one  of  the 
neighboring  slices.  After  the  procedure  is  completed,  the  welded  stack  geometry  is 

deleted  from  the  processor’s  memory  (together  with  the  corresponding  material  data). 

Near-field  data  in  the  compressed  matrix.  Finally,  in  order  to  generate  the  near- 
held  matrix  in  the  form  of  Eq.(3.7),  each  processor  p  splits  its  rectangular  (and  far-held 
subtracted)  matrix  block  into  separate  blocks,  A^p,p~a^  to  A^p'p+a\  with  columns 

Computation  of  the  mapping  coefficients  incurs  a  minor  cost,  hence  it  is  not  necessary  to  store  them  or 
transmit  them  between  the  processors. 
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indexed  by  unknowns  in  the  geometry  slices  Q^p  <J)  to  It  then  deletes  the  original 

block  A M. 


Communication  data  in  the  compressed  matrix.  The  compressed  matrix  compo¬ 
nents  constructed  as  described  above  have  to  be  complemented  with  information  on  how 
to  exchange  data  between  the  processors  in  matrix-vector  multiplications  executed  in  the 
iterative  solution  stage.  In  computing  both  near-  and  far-held  contributions  only  elements 
of  the  MoM  vectors,  such  as  x  of  Eq.(3.6),  are  being  exchanged.  Lists  of  such  elements 
( “communication  tables” )  are  constructed  in  the  preceding  steps  of  the  algorithm.  For  ex¬ 
ample,  each  processor  p  constructs  its  communication  tables  for  near-held  couplings  based 
on  the  structure  of  the  rectangular  near-held  matrix  block  A\p,a  1:  elements  of  the  vector 
blocks  x^p+T\  r/0,  that  have  to  be  received  from  other  processors  are  numbers  (3  of  those 
nonzero  columns  of  A^'^  which  reside  on  other  processors.  Similarly,  far-held  communica¬ 
tion  tables  are  built  on  the  basis  of  the  structure  of  the  sets  of  the  expansion  coefficients  V 
constructed  by  the  processor. 

3.2.3  Matrix-vector  multiplication 

Near-field  contributions.  With  the  data  organization  described  in  Section  3.1.4,  each 
processor  p  stores  a  row  (3.8)  of  near-held  matrix  blocks;  it  has  to  be  realized  that  the 
off-diagonal  blocks  A^p'p+T\  r  /  0,  may  only  contain  small  numbers  of  rows  and  columns. 
We  evaluate  then  the  p-th  block  of  the  near-held  matrix-vector  product  as 

a 

yNearfr)  ;=  A\p,a]  x  =  ^  A^’P+^  ,T(p+r)  (3.10) 

T=  —  (7 

(we  also  recall  here  that  a  has  been  dehned  as  the  maximum  number  of  neighboring  coupled 
slices;  the  actual  number  of  nonzero  blocks  in  the  set  (3.8)  may  be  smaller  than  2cr  +  1). 

In  order  to  compute  the  matrix  block  (3.10),  the  p-th  processor  has  to  receive  subsets  of 
the  vectors  x^p+T^  from  the  other  processors,  (p  +  r).  It  also  has  to  send  to  those  processors 
sub- vectors  of  its  own  input  vector  block  x^p\  These  operations  are  governed  by  the  near- 
held  communication  tables  constituting  a  part  of  the  compressed  matrix.  As  indicated  in 
Eq.(3.10),  they  are  executed  in  a  loop  through  the  index  r,  and  in  each  step  the  MPI 
function  MPI_Sendrecv  is  called  to  exchange  subsets  of  the  vector  x^p+T^\  those  elements 
that  are  being  sent  by  the  processor  p  have  been  previously  copied  to  and  stored  in  a  work 
buffer. 

Far-field  contributions.  This  part  of  the  algorithm  involves  the  operations  described 
by  Eqs.  (2.11)  to  (2.15). 

(i)  In  the  first  step,  each  processor  p  converts  its  MoM  vector  block  x ^  to  the  block  of 
the  Cartesian  representation  vector  X  associated  with  its  home  slice  C ^  of  the  Cartesian 
grid.  It  is  implemented  as 

x. ^  =  J2  E  V£+T'p)  xt+T)  for  a11  v  G  ciP)  •  (3.11) 

T=-<T  p<zQ(p+T) 
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The  (3  sum  runs  over  the  basis  functions  (unknowns)  in  the  geometry  slices  Q^p+T^  stored  on 
the  processor  p  itself  and  on  the  neighboring  processors  (p  +  r) .  The  expansion  coefficients 
y(p+T,p)  are  gtore(}  locally  on  the  processor  p,  but  in  a  sparse  representation  (Section  3.1.5), 
i.e.,  only  for  the  nodes  v  located  in  the  home  Cartesian  grid  slice  C^p\  As  in  the  near¬ 
field  matrix- vector  multiplication,  some  elements  of  the  vector  blocks  have  to  be  exchanged 
between  the  processors,  also  in  a  procedure  controlled  by  communication  lists  and  utilizing 
the  routine  MPI_sendrecv. 

(ii)  The  next  steps  involve  operations  (2.12),  (2.13),  and  (2.14),  on  the  Cartesian  vector 

A. 

The  forward  and  backward  FFTs  involve  global  data  transposition  (mentioned  in  the 
FFT  algorithm  description  in  Section  3.1.1),  i.e.,  global  communication  between  all  the 
processors.  Therefore,  its  implementation  is  crucial  for  the  parallel  performance  of  the 
solution  algorithm. 

In  our  solver  the  forward  FFT  (Eq.(2.12))  is  carried  out  with  the  help  of  FFTW  routines, 
but  with  a  communication  scheme  reducing  the  amount  of  the  transferred  data  (relative  to  a 
literal  application  of  the  parallel  FFTW  implementation  of  the  global  FFT).  Our  algorithm 
(described  in  Appendix  A)  carries  out  communication  by  directly  executing  the  MPI  vector 
all-to-all  communication  routine  MPI_Alltoallv. 

Next  (Eq.(2.13)),  each  processor  p  multiplies  its  segment  of  the  Fourier-transformed 
Cartesian  data  X  by  the  FFT  of  the  Green  function: 

y(p)  _  n<{p)  y(p) 

y'-nxnynz  ^nx  ny  nz  yynx  ny  nz  (3  12) 

for  0  <  nx  <  Nx  ,  0  <ny<  ,  0  <  nz  <  N z  . 

We  exhibited  here  explicitly  operations  on  the  indices  to  emphasize  the  fact  that  the  vec¬ 
tor  X  includes  padding  and  that  its  is  stored  in  the  transposed  representation,  and  thus 
distributed  across  the  processors  according  to  the  y-slices  (Eq.(3.4)). 

Finally,  the  backward  FFT  (Eq.(2.14))  is  executed  in  close  analogy  to  the  forward  trans¬ 
form,  with  the  reversed  order  of  its  operations. 

(iii)  The  last  operation  in  the  far-held  multiplication  is  the  conversion  of  the  held  data 
from  its  Cartesian  to  MoM  representations  (Eq.(2.15)).  It  is  implemented  in  parallel  as 

y  Far(p+r)  _  ^  V^rT'p">  for  t  =  —  <7,  ...  ,  a  and  for  all  a  £  Q ^  .  (3.13) 

ueC(p) 

As  in  Eq.(3.11),  the  expansion  coefficients  I/(p+t>p)  are  locally  stored  on  the  processor  p  in 
their  sparse  representation.  Also,  as  before,  the  procedure  involves  communication  between 
the  nearby  processors:  The  p-th  processor  evaluates  parts  of  the  vector  blocks  y(farip+ri, 
r  /  0,  and  sends  their  elements  to  the  processors  (p+r),  where  they  are  added  to  the  locally 
computed  elements.  At  the  same  time,  the  processor  p  receives  additional  contributions  to 
its  far-held  vector  block  from  the  processors  (p  +  r),  r  /  0. 

As  we  mentioned  in  Section  3.1.1,  we  the  far-held  matrix-vector  product,  Eq.(3.2),  as  a 
sequential  sum  of  contributions  of  the  four  components  of  the  Cartesian  equivalent  sources 
(the  pressure  and  the  three  components  of  its  gradient).  Therefore,  we  only  need  storage  for 
a  single-component  vector  X,  for  the  scalar  Green  function  G,  and  for  an  additional  work 
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buffer  used  in  the  (not  in-place)  communication  in  the  FFTs,  as  discussed  in  Appendix  A 
[...  not  really  discussed  now,  add  this];  the  size  of  the  buffer  is  one-half  of  the  size  of  X  or 
G. 

Computing  contributions  to  the  far-held  matrix-vector  product  sequentially  amounts  to 
transmitting,  in  the  parallel  FFT,  four  messages  each  of  the  quarter  size  of  the  full  message. 
It  might  seem  that  such  a  procedure  would  be  slower  than  transmitting  a  single  message  of 
the  same  total  size.  However,  in  our  computations  we  rather  observed  a  certain  speed-up; 
this  behavior  appears  to  be  due  to  the  low  latency  of  the  InfiniBand  interconnect  network 
and,  perhaps,  to  an  overall  positive  effect  of  a  significant  reduction  in  the  sizes  of  processes. 

4  Computational  complexity  estimates 

As  we  discussed  in  Ref.  [2] ,  the  complexity  of  the  serial  algorithm  algorithm  implementation, 
apart  of  matrix  fill,  is  dominated  by  the  FFTs  in  the  fast  matrix-vector  multiplication  in  the 
iterative  solution  stage.  As  follows  from  the  description  in  Section  3.2.2,  matrix  construction 
requires  a  minimal  amount  of  inter-processor  communication  and  is  expected  to  scale  with 
the  number  of  processors  practically  perfectly  -  a  feature  confirmed  by  results  of  many 
computations. 

The  near-field  part  of  the  matrix-vector  multiplication  in  the  iterative  solution  is,  for 
large  problems,  practically  negligible  in  comparison  the  far-field  part  of  the  algorithm.  The 
latter,  MoM  data  exchanges  between  nearby  slices  contribute  little  to  the  computational 
cost.  The  main  part  of  the  communication  overhead  is  associated,  as  we  discussed  in 
Section  3.2.3,  with  collective  (all-too-all)  communication  in  the  data  transposition  stage. 
However,  we  found  that  on  the  presently  available  parallel  cluster  systems  utilizing  fast 
interconnect  networks,  the  parallel  FFT  (i.e.,  in  fact,  the  collective  MPI  MPI_Alltoallv 
operation)  scales  quite  well:  although  the  speed-up  factor  is  not  proportional  to  P,  it 
scales,  typically,  as  P2//'3  to  P3/4.  As  a  result,  the  computational  cost  of  matrix- vector 
multiplication  is  quite  low.  For  instance,  in  one  of  our  examples  (Section  6)  it  amounts  to 
6.6  s  for  a  problem  with  N  ~  4,  700, 000  unknowns,  solved  on  128  processors. 

5  Storage  estimates 

In  view  of  the  effectiveness  of  parallelization  in  the  fast  matrix-vector  multiplication,  our 
main  concern  is  storage  per  processor,  which,  on  most  present  systems,  is  limited  to  mmax  = 
1.75  GB.  Therefore,  we  give  a  rather  detailed  discussion  of  the  storage  requirements  arising 
in  our  approach;  these  results  will  indicate  how  large  problems  can  be  solved  on  how  many 
processors,  given  the  existing  memory  restrictions. 

We  first  consider  a  rather  general  volumetric  problem,  and  express  the  required  storage 
in  terms  of  a  few  parameters  which  depend  on  the  object  geometry  and  its  discretization, 
and  have  to  be  estimated  for  specific  cases.  We  then  assume  a  uniform  and  perfectly 
balanced  distribution  of  the  data  and  obtain  simple  estimates  in  terms  of  the  geometry 
size,  its  discretization,  and  compression  parameters.  Finally,  we  discuss  limitations  on  the 
problem  size  imposed  by  the  approach  based  on  geometry  partition  into  slices  and,  possibly, 
on  the  use  of  the  interleaved  data  layout. 
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In  the  following  we  always  give  expressions  for  storage  on  a  single  processor,  say  p,  out 
of  P  >  1  processors. 


A  general  data  distribution.  We  assume  the  scatterer’s  bounding  box  is  oriented  ac¬ 
cording  to  the  conditions  (3.1),  and  that  the  box  sizes,  Fq,  i  =  x,y,z,  are  large  relative  to 
the  Cartesian  grid  spacing, 

h<.Bx<By<Bz.  (5.1) 

We  also  consider  discretization  with  fairly  uniformly  sized  and  regular  tetrahedra  of  average 
edge  length  a,  hence  the  average  volume 


Ma) 


6  a/2  ~  8.5 


(5.2) 


Further,  we  denote  by  p  the  ratio  of  the  grid  spacing  to  the  tetrahedron  size,  i.e.,  assume 


h  =  pa  . 


(5.3) 


We  notice  here  that  the  tetrahedron  volume  (5.2)  is  equal  the  Cartesian  cell  volume  /i3 
when 

p  =  p0  =  (6  v/2)~1/3  ~  0.490  ;  (5.4) 

incidentally,  this  value  of  p  appears  to  be  nearly  optimal  from  the  point  of  view  of  minimizing 
the  overall  computational  cost  of  matrix- vector  multiplication  [2]. 

Under  the  assumptions  made  above,  the  numbers  of  Cartesian  grid  points  assigned  to  a 
single  processor  p  are 


K  -  Bi  -  Bi 
1  h  pa 


i  =  x,  y,  z, 


K(p)  = 


P 


paP 


(5.5) 


It  follows  that  the  storage,  in  bytes,  for  a  single-component  Cartesian  vector,  such  as  G  or 
X,  is 

M\G\  =M\X]  =  B*B?^z  .  64 B  ,  (5.6) 

l  j  l  j  t]6  a6  Jr 

where  the  coefficient  8  •  8  =  64  is  due  to  the  factor  8  =  23  for  zero-padding  and  8  bytes  for 
a  single-precision  complex  number.  As  we  discussed  in  Section  3.2.3,  we  only  need  storage 
for  the  Green  function  G,  a  single  vector  X,  and  a  communication  buffer  B.  As  follows 
from  the  discussion  in  Appendix  A,  the  size  of  the  buffer  is  one-half  of  the  size  of  G  or  X. 
Hence,  the  total  storage  for  the  Cartesian  data  is 

~  B  B  B 

M[G,X,B]  =(1  +  1  +  \)M[X]  =  ^3Jpz-160B  .  (5.7) 

Storage  for  the  near-field  matrix  ANear  is  given  by 

M  [ANear]  =  JV(p)  N$  •  12  B  ,  (5.8) 
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where  N ^  is  the  number  of  unknowns  assigned  to  the  processor,  N ^  is  the  average  number 
of  nonzero  elements  per  row  of  the  matrix,  and  the  storage  of  12  B  per  matrix  element  arises 
from  the  sparse-row  matrix  representation  (a  complex  value  and  an  integer  offset). 

Clearly,  for  general  geometry,  the  number  of  unknowns  N ^  may  vary  in  a  wide  interval. 

(p) 

The  number  N ^  of  nonzero  elements  grows  proportionally  to  the  cube  of  the  near-field 
range;  its  typical  values  may  be  of  the  order  100  to  200. 

Finally,  we  give  an  estimate  for  the  far-held  part  of  the  compressed  stiffness  matrix, 
i.e.,  for  the  MoM-to-Cartesian  mapping  coefficients  V.  According  to  Eqs.  (2.7)  or  (3.2),  for 
each  relevant  tetrahedron  we  need  to  store  (1  +  3  +  3)  =  7  components  of  the  coefficients 
y(M),  V(D\  and  As  we  discussed  in  Section  3.1.5,  those  coefficients  have  to  be  stored 

for  all  tetrahedra  in  the  home  slice  of  the  processor,  and  for  some  nearby  tetrahedra  in 
the  adjacent  slices.  On  the  other  hand,  for  tetrahedra  near  the  slice  boundary,  we  store 
the  coefficients  not  for  all  nodes  of  the  expansion  boxes,  but  only  for  those  in  the  home 
Cartesian  slice.  In  fact,  as  we  find  below  (under  the  assumption  of  uniform  discretization) 
these  effects  tend  to  approximately  cancel,  with  the  result 

M[V\  ~JV(p)  -3024B  ,  (5.9) 

where  the  coefficient  7  •  (2  +  l)3  •  16  B  =  3024  B  is  obtained  for  M  =  2  and  for  complex 
double-precision  values  of  Vs. 


A  uniform  data  distribution.  We  consider  now  a  volumetric  rectangularly  shaped  ob¬ 
ject  completely  filling  its  bounding  box  of  sizes  (5.1),  and  uniformly  discretized  with  tetra¬ 
hedra  of  edge  length  a. 

Under  these  simplifying  assumptions,  the  number  of  tetrahedra  (unknowns)  per  proces¬ 
sor  is  given  by  the  ratio  of  the  slice  volume  to  the  tetrahedron  volume  (5.2),  i.e., 


N(p) 


Bx  By  Bz 

r] o  a3  P 


~  8.5 


Bx  By  Bz 

a3  P 


(5.10) 


Similarly,  the  number  of  nonzero  elements  per  row  of  the  near-field  matrix  is  given  by  the 
ratio  of  the  volume  of  a  ball  of  radius  i?N  and  the  tetrahedron  volume, 


iy(p) 

ivN 


4vr  R3 j 
3  rfa  a3 


~  35.5 


Ri 


'N 


(5.11) 


We  note  that  the  typical  compression  parameters  (r/  =  rj0  and  i?N  =  3 h)  result  in  = 
36  7r  ~  113. 

Eqs.  (5.10)  and  (5.10)  yield  now  an  estimate  of  the  near-field  matrix  storage  (5.8)  as 


T  Nparl  BxB  B  Rfi  B  B  B,R\ t 

M  ANear  =  32  7T  l  *  N  •  12  B  ~  yr  ~  N  •  1206  B  .  5.12 

L  J  a6  P  a6  P  K 

For  the  considered  uniform  geometry  discretization  we  can  also  obtain  a  more  precise 
estimate  of  the  storage  for  the  mapping  coefficients  V ,  i.e.,  for  the  number  of  grid  nodes, 
associated  with  the  relevant  tetrahedra,  for  which  the  mapping  coefficients  have  to  be  stored 
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(as  discussed  in  Section  3.1.5).  To  facilitate  the  reasoning,  we  first  assume  that  the  centroids 
of  the  tetrahedra  coincide  with  the  grid  nodes  (or  the  grid  cell  centroids),  and  then  multiply 
the  result  by  the  factor  (?//?7o)3,  which  accounts  for  the  difference  in  the  distribution  density 
of  the  tetrahedra  and  grid  nodes. 

We  assume  in  the  following  M  >  1 ,  and  first  consider  “thick  slices” ,  by  which  we  mean 

(p) 

Kz  >  M  +  1.  In  this  case  there  exist  expansion  boxes  entirely  contained  in  the  Cartesian 
home  grid  slice  C^p\  By  counting  the  numbers  of  grid  nodes  u  located  in  the  home  grid 
slice  for  expansion  boxes  in  the  home  slice  and  intersecting  its  boundaries,  and  multiplying 
the  result  by  the  volume  factor,  we  obtain  the 

#u  =  (— \  Kx  Ky  [l4p)  (. M  +  1)  -  M  +  2]  (M  +  l)2  for  K {p)  >M+  1  >  2  .  (5.13) 
\VoJ 

(p) 

For  “thin  slices”,  Kz  <  M,  a  similar  reasoning  yields  simply 

#u  =  I\x  Ky  /dp)  (AT  +  l)3  for  K<f>  <  M  .  (5.14) 

\VoJ 

In  the  usual  case  M  =  2,  both  of  the  above  expression  reduce,  in  view  of  Eqs.  (5.5)  and 
(5.10),  to  the  simple  expectation 

#u  =  Nip)  (M  +  l)3  for  M  =  2  ,  (5.15) 


which  immediately  implies  Eq.(5.9). 

The  above  storage  estimates  for  the  mapping  coefficients  V  apply  to  the  “sparse”  V  rep¬ 
resentation,  in  which,  for  every  tetrahedron,  the  coefficients  are  stored  only  for  the  relevant 
Cartesian  nodes  (those  in  the  home  grid  slice).  Were  we  to  use  the  “full”  representation, 
i.e. ,  store,  for  every  contributing  tetrahedron  a,  the  mapping  coefficients  for  all  the  (M  +  1)3 
nodes  of  the  expansion  box  B(a),  the  storage  would  have  been,  instead  of  Eq.(5.9), 

X  [r,  J  -  N<r)  (l  +  H) )  '  302“  B  ■  (5.16) 

\  J\Z  / 

The  savings  due  to  the  sparse  storage  are  thus  significant,  especially  for  thin  slices:  if 

(p) 

Kz  '  <  M,  then  the  ratio  of  the  full  to  sparse  storage  size  is  at  least  3, 

2  M 

l+^S3.  (5.17) 

(P) 

and  may  reach  5  in  the  extreme  case  of  K'/  =  1,  and  for  M  =  2. 

To  summarize,  under  the  simplifying  assumptions  made  above,  and  with  the  resulting 
expression  (5.10)  for  the  number  of  unknowns,  the  total  storage  can  be  represented  as 
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(5.18) 


where  N  is  the  total  number  of  unknowns.  As  follows  from  Eqs.  (5.12),  (5.9),  and  (5.7), 
the  amount  of  storage  per  unknown  is  independent  of  the  problem  size  and  given  by 
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where  the  three  terms  are  due,  respectively,  to  the  near-field,  the  mapping  coefficients  V, 
and  the  Cartesian  vector  storage.  The  reason  we  expressed  the  estimate  (5.19)  in  terms  of 
R^/h  rather  than  i?N/a,  is  that  the  accuracy  of  the  matrix  compression  (2.5)  depends,  at 
least  in  the  typical  range  of  parameters,  mostly  on  the  first  ratio  [2].  Therefore,  storage 
can  be  minimized  by  varying  (although  in  a  limited  range)  the  parameter  77  while  keeping 
R^/h,  and  hence  the  accuracy,  fixed.  Without  entering  into  the  details,  we  can  state  that, 
depending  on  the  material  parameters  and  thus  the  magnitude  of  various  terms  in  the 
integral  equation  kernel,  the  required  value  of  RN/h  may  vary  from  about  3  to  about  6, 
while  i]  is  limited  to  i]  >0.6  ?70 -  Correspondingly,  the  storage  m0  may  vary  in  a  wide  range, 
from,  say,  5kB  to  20  kB. 

Finally,  we  comment  here  on  the  potential  disadvantages  of  our  “interleaved”  Cartesian 
storage  scheme,  mentioned  in  the  Introduction.  We  note  that,  for  the  range  of  parameters 
considered  above,  the  entire  Cartesian-vector  storage  (the  last  term  in  Eq.(5.19))  takes  at 
most  about  one-quarter  of  the  V  storage,  and  even  less  of  the  total  storage,  especially  for 
larger  near-field  ranges  (about  10%  for  R^/h  =  6).  Out  of  this  fraction,  2/5  constitutes 
additional  storage  for  padding  (the  drawback  (a))  and  1/5  the  storage  for  the  communication 
buffer  (the  drawback  (b)).  Thus,  we  can  conclude  that  the  additional  storage  cost  in  the 
interleaved  data  layout  is  minor  and  should  be  easily  compensated  by  the  advantages  of  the 
scheme. 


Limitations  on  the  problem  size.  The  maximum  size  of  treatable  problems  is,  obvi¬ 
ously,  limited  by  the  condition  A4tot  <  rnmax,  or,  in  view  of  Eq.(5.18),  by 

/V  Tfl 

4-  <  ^222*  (5.20) 

P  m0 

a  value  in  the  range  from  about  90, 000  to  370, 000. 

At  the  same  time,  however,  a  scheme  based  on  geometry  partitioning  into  slices  does 
not  allow  us  to  use  more  processors  than  the  total  number  of  Cartesian  grid  points  in  the 
z  direction,  Kz.  The  severity  of  this  restriction  depends,  of  course,  on  the  geometry  of 
the  object.  In  the  simplest  case  of  a  cubic  shape,  with  Bx  =  By  =  Bz  =  B,  we  can 
use  Eqs.  (5.10),  (5.5),  and  (5.3)  to  relate  Kz  to  the  number  of  unknowns,  and  obtain  the 
condition 

P  <  K  =  %  =  ^  AT1/3  .  (5.21) 

h  Tj 

Eqs.  (5.20)  and  (5.21)  provide  now  lower  and  upper  bounds  on  the  number  of  processors, 

— N  <  P  <  —  N1/3  ,  (5.22) 

^®max  V 

which  impose  a  limit  on  the  number  of  unknowns, 

N<(!h  Hhs**}  '  .  (5.23) 

V  V  mo  ) 

Thus,  for  m0  =  20  kB  and  rj  =  0.6  ?70 ,  and,  under  the  idealized  conditions  we  assumed, 
N  <  60,000,000. 
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6  Examples  of  large-scale  computations 

6.1  The  geometrical  models  and  material  properties 

We  present  below  results  of  example  computations  involving 


1.  a  realistically-shaped  model  of  a  human  head; 

2.  the  same  head  model  and  a  model  of  helmet,  with  the  in-between  space  filled  by  air; 

3.  the  model  of  a  head  and  a  helmet,  with  the  in-between  space  filled  by  a  low-density 
material  (cork). 

We  used  two  tetrahedral  discretizations  of  the  models,  differing  by  about  a  factor  2  in  the 
linear  sizes  of  the  tetrahedra;  the  average  edge  lengths  were  about  6.3  mm  and  3.3  mm.  The 
models  of  the  head  and  the  helmet,  with  the  finer  discretization,  are  shown  in  Fig.  2,  and 
the  parameters  of  the  geometries  are  listed  in  Table  1.  The  numbers  given  there  indicate 
a  fairly  uniform  discretization:  in  the  worst  case  the  minimum  edge  length  is  about  1  / 5-th 
of  the  average  value,  and  the  minimum  angle  of  a  facet  is  4.3°.  Volumetric  discretization 
(tetrahedronization)  of  surface  models  was  carried  out  by  means  of  the  program  tetgen, 
version  1.4.1,  available  in  the  public  domain. 


Figure  2:  Model  of  the  external  head  surface  and  the  helmet.  The  arrow  indicates  the 
direction  of  the  incident  wave. 
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Table  1:  Parameters  of  the  geometries  used  in  the  computations:  the  number  of  vertices 
v,  the  number  of  edges  e,  the  number  of  facets  /,  the  number  of  tetrahedra  t,  average 
edge  length  £av ,  minimum  edge  length  ^min  (both  in  meters) ,  minimum  facet  angle  amin  (in 
degrees) 


geometry 

V 

e 

/ 

t 

f 

^av 

f  . 

^min 

rv  • 

head 

65,797 

471,186 

803,108 

397,718 

0.00626 

0.00233 

16.8 

head  +  helmet 

82,708 

556,731 

923,469 

449,444 

0.00619 

0.00221 

16.8 

head  +  cork  +  helmet 

112,510 

801,487 

1,364,352 

675,374 

0.00627 

0.00221 

16.8 

head 

442,065 

3,214,434 

5,512,045 

2,739,675 

0.00329 

0.00112 

13.3 

head  +  helmet 

516,465 

3,593,613 

6,047,750 

2,970,600 

0.00331 

0.00059 

4.3 

head  +  cork  +  helmet 

754,796 

5,467,095 

9,366,177 

4,653,875 

0.00331 

0.00059 

4.3 

We  assumed  the  head  model  filled  with  a  homogeneous  material,5  with  mechanical 
properties  of  bone  (which  would  effectively  constitute  the  the  outer  layer  of  the  head) .  The 
acoustic  parameters  of  the  materials  are  listed  in  Table  2. 

The  models  were  placed  in  the  field  of  a  harmonic  acoustic  wave  of  frequency  5  kHz  (and 
the  wavelength  in  air  A  =  6.8 cm),  incident  laterally  on  the  right  ear  of  the  head  model,  as 
indicated  in  Fig.  2. 

The  computations  were  carried  out  on  a  Linux  cluster  with  the  InfiniBand  interconnect 
network,  on  32,  108,  and  128  processors. 


Table  2:  Acoustic  properties  of  materials  used  in  the  simulations 


material 

p/po 

n2 

bone 

1777.0 

0.1524 

brain 

835.4 

0.0564 

cork 

150.0 

0.65 

steel 

6667.0 

0.0035 

In  the  following  Subsections  we  discuss  some  aspects  of  the  computation,  such  as  the 
accuracy  of  the  far- field  matrix  compression,  physical  properties  of  the  solutions,  and  scaling 
of  the  code  performance. 

5We  stress,  however,  that  the  presence  of  inhomogeneities  would  not  increase  the  computational  cost  of 
the  solution. 
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6.2  Analysis  of  the  accuracy  of  the  far-held  expansion  in  the  problem  of 
the  human  head  helmet  with  cork  filling 

The  accuracy  of  the  compressed  representation  of  the  far-held  depends  on  two  compression 
parameters:  the  Cartesian  grid  sampling  ratio  y  :=  A /h  and  the  near-held  range  (in  the 
solution  of  the  surface  problem)  measured  in  units  of  the  Cartesian  grid  spacing,  i.e. ,  the 
parameter  R^/h.  We  present  here  the  comparison  of  solutions  obtained  for  a  hxed  value 
of  the  sampling  ratio,  y  :=  X/h  (i.e.,  a  hxed  Cartesian  grid  spacing  h),  and  several  values 
of  the  range  parameter,  i?N//i  =  6,8, 10, 12.  In  all  cases  the  near-held  range  in  the  volume 
problem  solution  is  one-half  that  in  the  surface  problem. 

As  we  mentioned  following  Eq.(5.19),  the  far-held  accuracy  is  controlled  primarily  by 
the  value  of  the  ratio  i?N//i,  and  depends  only  weakly  on  the  Cartesian  grid  spacing  h  (in 
a  reasonable  range  of  values  h<a,  where  a  is  the  average  tetrahedron  size).  Therefore,  by 
keeping  h  constant  and  varying  R^/h,  we  are  essentially  testing  the  sensitivity  of  solutions 
to  the  far-held  compression  accuracy. 

As  we  discussed  in  Sec.  5,  the  computational  cost  of  the  solution  (more  precisely,  the 
matrix  fill)  and  the  required  storage  strongly  depend  on  the  range  parameter  f?N.  In 
particular,  the  size  of  the  near-held  component  of  the  stiffness  matrix  is  approximately 
proportional  to  R j(,  (Eqs.  (5.11),  (5.12),  and  (5.19)). 

The  analysis  was  carried  out  for  the  largest  problem  we  considered,  i.e.,  the  model  of  a 
human  head  with  a  helmet,  with  the  space  between  the  two  filled  with  cork.  As  indicated 
in  Table  1,  the  system  was  discretized  with  N  ~  4,  700,  000  tetrahedra.  The  solutions  were 
computed  on  P  =  128  processors. 

Fig.  3  shows  the  distribution  of  the  far-held  errors  £aa  (Eq.(2.9))  as  a  function  of  the 
distance  R  =  Rap  =  |cQ  —  c^|  between  the  centroids  of  the  basis  functions  <pa  and  (f>p 
(Eq.(2.8)).  More  precisely,  the  range  of  distances  R  is  divided  into  nearly  200  bins,  and  in 
each  bin  the  average  error  is  computed. 

Fig.  3  indicates  that,  in  order  to  achieve  a  given  accuracy,  the  range  R  must  be  signifi¬ 
cantly  (about  a  factor  of  2)  larger  in  the  surface  problem  than  in  the  volume  problem.  This 
behavior  is  expected  on  theoretical  grounds,  since  the  matrix  elements  in  the  volume  prob¬ 
lem  are  dominated  by  monopole-monopole  interactions,  while  those  in  the  surface  problem 
involve  only  monopole-dipole  couplings.  Therefore,  in  the  latter  case,  the  dipole  held  has 
to  be  approximated  by  a  set  of  monopole  sources  located  in  the  “expansion  box”  of  the 
considered  basis  function,  and  such  an  approximation  cannot  be  valid  in  the  immediate 
vicinity  of  the  box. 

As  we  stated  above,  we  compare  results  of  computations  for  near-held  range  parameters 
RN/h  =  6,8,10,12  (referring  to  the  surface  problem;  in  the  volume  problem  we  assume 
ranges  equal  one-half  of  these  values).  For  those  choices,  the  error  values  at  the  maximum 
distance  drop  from  about  3%  to  below  1  %,  in  both  surface  and  volume  problems. 
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Raft/  h 

Figure  3:  Averaged  (over  distance  bins)  relative  far-field  errors  in  the  surface  (S)  and  volume 
(V)  problems  for  the  range  parameter  R^/h  =  12,  plotted  as  functions  of  the  distance 
between  basis  function  supports,  measured  in  units  of  the  Cartesian  grid  spacing  h.  The 
plots  represent  output  of  the  processor  p  =  59  with  the  assigned  ~  60, 000  unknowns. 

As  the  first  check  on  the  effects  of  the  near- field  range  parameter  we  show  in  Fig.  4 
plots  of  convergence  of  the  iterative  solutions.  It  is  evident  that  convergence  in  the  surface 
problem  is  practically  independent  of  i?N.  In  the  volume  problem  the  solution  for  R^/h  =  6 
converges  differently  than  in  the  remaining  cases,  but  for  R^/h  >  8  the  curves  are  nearly 
indistinguishable. 
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Figure  4:  Convergence  (the  relative  residual  norm)  of  iterative  solutions  in  the  surface  (S) 
and  volume  (V)  problems  for  the  indicated  values  of  the  near-field  range  parameter  R^/h. 


The  iterative  convergence  results  already  suggest  that  even  the  smallest  near-field  range 
and  the  resulting  far-held  error  on  the  level  of  ~  3%  are  sufficient  for  obtaining  accurate 
solutions.  This  assertion  is  verified  by  the  following  plots  of  pressure  distributions. 

Fig.  5  shows  pressure  distributions  on  the  exterior  surface  of  the  model.  These  results, 
obtained  by  solving  only  the  surface  problem,  are  visually  practically  indistinguishable. 


(a)  (b)  (c)  (d) 


Figure  5:  Distributions  of  the  real  part  of  pressure  on  the  external  surface  of  the  model 
of  the  human  head,  the  helmet,  and  the  cork  Filing,  with  N  ~  4,700,000  unknowns.  The 
distributions,  plotted  in  linear  scale,  were  obtained  for  the  same  values  of  the  near-held 
range  as  in  Fig.  4,  i.e. ,  (a)  RN/h  =  6,  (b)  RN/h  =  8,  (c)  R^/h  =  10,  (d)  R^/h  =  12. 

The  following  Figures,  6  and  7,  show  pressure  distributions  inside  the  model.  These 
solutions  are  obtained  from  the  volumetric  integral  equations  in  the  second  step  of  the 
solution  procedure,  which  uses  the  surface  problem  solution  as  input. 
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In  order  to  exhibit  the  oscillatory  nature  of  the  solution  we  plotted  in  Fig.  7  the  real  part 
of  the  pressure.  For  additional  emphasis,  that  distribution  is  displayed  in  the  logarithmic 
scale  (modified  near  zero  values).  The  pattern  of  oscillation  on  the  object  surface  matches 
that  shown  in  Fig.  5. 


16.0 


8.0 


o.o 


Figure  6:  Distributions  of  the  absolute  value  of  pressure  in  the  coronal  section  of  the  model 
of  the  human  head,  the  helmet,  and  the  cork  filling,  with  N  ~  4,700,000  unknowns.  The 
distributions,  plotted  in  linear  scale,  were  computed  for  the  same  values  of  the  near-field 
range  as  in  the  previous  Figures,  i.e. ,  (a)  R^/h  =  6,  (b)  =  8,  (c)  R-^/h  =  10, 

(d)  R^/h  =  12. 
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Figure  7:  Distributions  of  the  real  part  of  pressure  in  the  coronal  section  of  the  model  of 
the  human  head,  the  helmet,  and  the  cork  filling,  with  N  ~  4,  700,  000  unknowns.  The 
distributions,  plotted  in  logarithmic  scale,  were  obtained  for  for  the  same  values  of  the 
near-field  range  before,  i.e.,  (a)  R^/h  =  6,  (b)  i?N//i  =  8,  (c)  R^/h  =  10,  (d)  R^/h  =  12. 

To  summarize,  Figs.  3  to  7  indicate  that  the  near-field  range  R^/h  =  8  (or  even  R^/h  = 
6),  resulting  in  far-held  expansion  accuracy  on  the  level  of  2%  (or  3%),  yields  results 
practically  indistinguishable  from  those  for  a  much  large  range,  R^/h  =  12. 
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6.3  Pressure  distributions  in  a  human  head  model:  dependence  on  the 
discretization 

We  compare  now  pressure  distributions  in  the  head  (without  the  helmet),  computed  for  two 
discretizations.  In  the  more  coarsely  discretized  model  (a)  the  number  of  tetrahedra  (and 
unknowns)  is  IV  ~  400,000;  the  finer  discretization  (b)  yields  TV  ~  2,  700,000  (Table  1). 

The  computations  have  been  carried  out  with  the  Cartesian  grid  sampling  rate  x  =  80 
and  with  the  near-field  range  i?N//i  =  6.  As  we  have  discussed  in  Subsection  6.2,  these 
compression  parameters  are  entirely  sufficient  for  obtaining  accurate  solutions.  The  total 
solution  times  were  (a)  about  30  minutes,  with  2.5  s  per  iteration,  on  P  =  32  processors, 
and  (b)  about  50  minutes  with  4.6  s  per  iteration,  on  P  =  108  processors. 

The  results,  shown  in  Fig.  8,  demonstrate  a  good  agreement  between  the  pressure  dis¬ 
tributions  obtained  for  the  two  discretizations. 
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Figure  8:  Distributions  of  the  real  part  of  the  pressure,  plotted  in  the  logarithmic  scale,  on 
several  sections  in  the  axial  plane  of  the  human  head  model  with  (a)  a  coarser  discretization 
( N  ~  400,  000  tetrahedra),  and  (b)  a  finer  discretization  (N  —  2,  700, 000  tetrahedra).  The 
models  are  subject  to  an  acoustic  wave  of  unit  pressure  amplitude  and  frequency  5  kHz, 
incident  horizontally  on  the  right  ear.  In  (b)  the  outline  of  the  head  outer  surface  is 
superimposed  on  the  pressure  distributions. 

6.4  Pressure  distributions  in  a  human  head  model  in  the  presence  of  an 
unattached  helmet 

Next,  we  show  results  for  the  human  head  model  with  the  finer  discretization,  suspended 
inside  a  helmet.  The  problem  involves  N  —  3,000,000  (Table  1)  and  has  been  solved  on 
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P  =  108  processors,  with  the  same  matrix  compression  parameters  as  before  (x  =  80, 
RN/h  =  6). 

The  results  (Fig.  9)  show  the  pressure  distribution  inside  the  head  roughly  similar,  and 
of  similar  magnitude,  as  for  the  isolated  head  (Fig.  8(b)),  although  the  details  are  different. 
One  can  thus  infer  that  the  presence  of  a  helmet  -  not  attached  in  any  way  to  the  head 
-  has  only  a  moderate  effect  on  the  pressure  values  inside  the  head  model.  However,  a 
configuration  with  a  helmet  “floating”  in  the  air  can  hardly  be  considered  a  physically 
realistic  model. 
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Figure  9:  Distribution  of  the  real  part  of  the  pressure,  plotted  in  the  logarithmic  scale,  on 
several  sections  in  the  axial  plane  of  the  human  head  model  with  an  unattached  helmet, 
and  the  total  number  of  unknowns  N  ~  3, 000,  000.  The  frequency  and  the  wave  incidence 
direction  are  as  in  previous  Figures. 

6.5  Pressure  distributions  in  a  human  head  model  in  the  presence  and 
absence  of  a  helmet  with  a  cork  lining 

As  the  final  example,  we  compare  the  pressure  distribution  in  the  head  model  obtained 
in  the  previous  calculation  (for  the  finer  discretization)  with  the  distribution  for  the  same 
model  with  an  addition  of  a  steel  helmet  and  the  in-between  space  filled  with  cork.  The 
latter  geometry  is  discretized  with  IV  ~  4,  700,  000  tetrahedra,  also  with  the  tetrahedron 
size  of  about  3  mm  (Table  1).  The  solution  for  the  head  and  helmet  system  was  carried  out 
on  128  processors  in  the  total  time  of  about  2  hours,  and  6.6  s  per  iteration.  The  overall 
time  increase  in  that  case  is  mostly  due  to  the  increased  number  of  iterations. 

We  have  already  described  a  set  of  computations  for  the  head,  cork,  and  helmet  system 
in  Subsection  6.2.  Here  we  comment  on  the  physical  features  of  the  problem  and  its  solution, 
and  their  relation  to  the  problem  of  the  isolated  head. 
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The  results  of  the  computations  are  visualized,  as  distributions  of  the  absolute  value  of 
the  pressure  in  the  coronal  plane  of  the  models,  in  Fig.  10.  Fig.  10(ii)  is  identical  to  the 
previous  Fig.  6(b).  The  real  part  of  the  pressure  has  been  shown  before  in  Fig.  7(b). 

The  solutions  show  a  nontrivial  behavior  and  exhibit  physical  phenomena  which  may 
be  relevant  in  the  design  of  protective  devices: 

In  the  case  of  the  head,  Fig.  10(i),  the  pressure  is  maximal  at  the  entrance  to  the  ear 
canal,  and  it  is  smoothly  distributed  inside  the  head:  while  the  wavelength  in  the  air  is 
about  6.8  cm,  it  is  6. 8cm /n  =  6. 8cm /\/0.1524  =  17.4  cm  in  the  head  material  (it  would 
have  been  even  longer  in  soft  tissues,  e.g.,  in  the  brain).  In  fact,  the  solution  is  suggestive 
of  a  resonance-type  (P-wave)  behavior:  the  pressure  changes  sign  along  the  approximately 
vertical  line  seen  in  the  Figure  and,  more  distinctly,  in  the  spatial  distribution  shown  in 
Figs.  8. 

The  solution  for  the  head  and  helmet  system,  Fig.  10(ii),  is  quite  different.  It  exhibits 
a  distinct  oscillatory  behavior  along  the  surface  of  the  helmet  and  in  the  region  filled  by 
cork.  This  region  appears  to  have  properties  of  a  “waveguide” :  because  of  the  cork  density 
being  significantly  lower  than  that  of  the  surrounding  materials  (the  helmet  and  the  head), 
and  the  resulting  impedance  mismatch  at  the  boundaries,  the  wave  tends  to  be  trapped  in 
that  region.  Since  the  refractive  index  of  cork  is  not  much  different  from  that  of  air,  wave 
oscillations  are  relatively  rapid.  We  stress,  however,  that  the  physical  picture  suggested  by 
Fig.  10(ii)  would  change  if  we  considered  a  dissipative  (attenuating)  filling  material,  e.g.,  a 
strongly  damping  porous  material  characterized  by  a  complex  refractive  index. 

We  note,  finally,  that,  for  the  particular  frequency  considered  here,  the  presence  of  the 
helmet  with  a  cork  lining  completely  changes  the  pressure  distribution  inside  the  head, 
but  does  not  reduce  its  maximum  value  (we  note  that  the  data  Figs.  10(i)  and  10(ii)  are 
rendered  in  different  scales).  This  feature,  however,  is  likely  to  depend  on  properties  of  the 
lining  (filling)  material. 
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(i)  (ii) 


Figure  10:  Pressure  distributions  in  the  coronal  plane  for  (i)  the  human  head  model  and  (ii) 
the  system  consisting  of  the  human  head  and  a  steel  helmet  models,  with  the  in-between 
space  filled  by  cork.  The  numbers  of  unknowns  in  these  two  problems  are  N  ~  2,700,00 
and  N  —  4,700,00.  The  maximum  pressure  values  are  about  4  in  (i)  and  15  in  (ii). 

6.6  Scaling  in  parallel  computations 

To  assess  scaling  of  our  solver  (with  the  problem  size  and  the  number  of  processors)  we 
summarize  the  results  of  a  set  of  computations  for  the  models  described  at  the  beginning 
of  this  Section.  We  list  in  Table  3  the  time  of  matrix- vector  multiplication6  as  a  function 
of  (a)  the  number  of  unknowns  IV;  (b)  the  Cartesian  grid  sizes  Kx,  I\y,  and  Kz  (always 
without  padding);  (c)  the  number  of  processors  P;  (d)  the  number  of  Cartesian  grid  nodes 
per  processor,  Kc/P  :=  Kx  I\  y  Kz/P ;  and  (e)  the  near-field  range  in  units  of  the  Cartesian 
grid  spacing,  RN/h  (which  affects  the  size  of  the  near-field  part  of  the  matrix). 

®More  precisely,  we  list  the  time  per  iteration,  which  includes  additional  computations  of  vector  dot- 
products  and  of  linear  combinations  of  vectors. 
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Table  3:  Matrix  compression  parameters  used  in  the  computations 


N 

Kx 

X 

Ky 

X 

Kz 

P 

KJP 

R^/h 

t  (s) 

397,718 

160 

X 

192 

X 

180 

32 

172,800 

3.2 

2.1 

397,718 

180 

X 

216 

X 

180 

32 

218,700 

6.0 

2.5 

449,444 

180 

X 

200 

X 

192 

32 

216,000 

3.2 

2.6 

449,444 

192 

X 

216 

X 

192 

32 

248,832 

6.0 

3.5 

675,374 

180 

X 

200 

X 

192 

32 

216,000 

4.8 

2.4 

2,739,675 

300 

X 

384 

X 

324 

108 

345,600 

6.0 

4.6 

2,970,600 

375 

X 

405 

X 

384 

128 

455,625 

6.0 

7.0 

4,653,875 

375 

X 

405 

X 

384 

128 

455,625 

6.0 

6.9 

4,653,875 

375 

X 

405 

X 

384 

128 

455,625 

12.0 

7.1 

Even  though  the  number  of  unknowns  and  the  number  of  processors  vary  by  the  factors 
of  12  and  4,  the  Table  shows  that  the  matrix-vector  multiplication  time  depends  mainly 
on  a  single  parameter:  the  number  of  Cartesian  grid  nodes  per  processor.  This  behavior  is 
expected  when  matrix-vector  multiplication  complexity  is  dominated  by  the  Fourier  trans¬ 
form,  which  also  approximately  scales  with  the  number  of  processors.  These  assumptions 
imply 

t(KC »  p)  -  7  -p  Kc  loS  2 (Kc)  (lo§  2^T  +  log  2P)  (6-1) 

with  some  fixed  coefficient  7  (depending  also  on  the  hardware,  particularly  on  the  commu¬ 
nication  network). 

Fig.  11,  in  which  we  plot  the  matrix-vector  multiplication  times  for  the  problems  of 
Table  3,  shows  that  the  actual  performance  of  the  solver  is  close  to  the  behavior  predicted 
by  Eq.(6.1). 
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scaling  of  matrix- vector  multiplication  time 


KC/P 


Figure  11:  Matrix-vector  multiplication  time  as  a  function  of  the  number  of  Cartesian  grid 
nodes  per  processor,  Kc/P,  for  the  computations  listed  in  Table  3.  The  curves  are  defined 
by  Eq.(6.1)  with  7  =  0.6  •  10-6  s,  and  plotted  for  P  =  32  (the  lower  curve)  and  P  =  128 
(the  upper  curve).  Numbers  of  processors,  P,  are  shown  next  to  the  corresponding  groups 
of  points. 
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Appendices 


A  Handling  of  zero  padding  and  communication  in  parallel 
implementation  of  FFTs 

We  discuss  here  some  aspects  some  aspects,  essential  for  the  code  performance,  of  parallel 
implementation  of  the  FFTs. 

In  Fig.  12  we  illustrate  the  operations  carried  out  on  the  Cartesian  vector  X  in  the  case 
of  P  =  3  processors.  In  the  problem  shown  in  the  Figure  the  Cartesian  vectors  are  specified 
in  terms  of  their  sizes  Ki  (without  padding)  and  Ni  =  2  Ki  (with  padding);  in  the  Figure 
we  do  not  display  the  x-dimensions,  and  leave  them  unspecified. 


Z 


(a) 


y 


(b) 


Figure  12:  FFT  in  the  (x,  y )  directions  of  the  distributed  data  (a),  followed  by  transposition 
and  FFT  in  the  z  direction  (b).  The  FFTs  extend  the  range  of  data  from  the  original 
storage  without  padding  (grey  background)  to  the  storage  with  padding  (grey  +  cross- 
hatched  areas).  Data  owned  by  the  individual  processors  are  shown  in  boxes  labeled  with 
the  processors  numbers  p.  The  zig-zag  lines  indicate  the  data  storage  order.  Further  details 
are  described  in  the  text. 


Our  parallel  FFT  algorithm  implementation,  while  it  uses  the  FFTW  routines,  is  sig¬ 
nificantly  more  efficient  than  the  original,  totally  “encapsulated”  three-dimensional  FFT 
provided  in  the  FFTW  package.  In  fact,  the  “encapsulated”  FFT,  with  the  transforma¬ 
tion  plan  generated  by  the  routine  fftw_plan  f ftw_mpi_plan_dft_3d  is  practically  not 
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usable  in  our  solver.  The  principal  reason  is  a  conflict  between  the  way  this  routine  handles 
padding,  and  an  economical  distribution  of  data  in  the  entire  matrix-vector  multiplication 
scheme: 

FFT  padding  in  the  z  direction  has  to  occupy  a  contiguous  range  of  nz  indices;  hence, 
if  z-slices  of  the  Cartesian  grid  are  distributed  across  the  processors,  about  one-half  of 
them  will  own  exclusively  the  padding  space.  This  configuration  is  more-or-less  adequate 
if  only  FFTs  are  being  computed,  although  even  in  this  case  handling  of  zero  padding 
is  not  efficient.  However,  the  other  operations  in  the  matrix- vector  multiplication  (near¬ 
field  computation  and  conversion  between  the  MoM  and  Cartesian  representations)  require 
access  to  the  MoM  vector  x  and  to  the  multipole  expansion  coefficients  V,  and  question 
arises  how  to  to  distribute  these  data.  There  are  at  least  two  possibilities: 

1.  Distribute  slices  of  the  x  and  V  data  only  across  the  processors  storing  the  “physical” 
segment  of  the  Cartesian  grid.  In  this  case  near-field  and  MoM-Cartesian  conversion 
operations  can  be  performed  locally,  but  only  one-half  of  the  processors  participates. 

2.  Distribute  slices  of  the  x  and  V  data  across  all  processors.  In  this  arrangement  all 
processors  are  involved  in  the  near-field  computation  and  MoM-Cartesian  conversion, 
but  the  resulting  Cartesian  vector  data  have  to  be  redistributed  before  and  after 
performing  the  FFTs. 

Evidently,  none  of  them  is  computationally  efficient. 

Another  potential  source  of  inefficiency  of  the  encapsulated  three-dimensional  FFTW 
transform  is  the  redundant  (in  our  case)  “back-transposition”,  which  the  FFTW  routine 
performs  in  order  to  restore  the  original  layout  of  the  data.  This  operation  can  be  eliminated 
by  creating  the  forward  and  backward  FFT  plans  with  the  flags  FFTW_MPI_TRANSPOSED_OUT 
and 

FFTW_MPI_TRANSPOSED_IN,  which  suppress  the  “back-transposition”,  and  cause  the  Fourier- 
transformed  data  to  be  left  in  the  transposed  storage  form.  Nevertheless,  the  previously 
listed  difficulties  remain. 

The  procedure  we  implemented  in  our  code  avoids  the  difficulties  entirely:  both  the  FFT 
and  the  {x,  V}  data  are  distributed  across  all  processors,  all  processors  participate  approx¬ 
imately  equally  in  the  near-  and  far-held  operations,  as  well  as  in  the  FFT  computations, 
and,  besides,  by  the  appropriate  data  arrangement,  we  reduce  communication  costs  by  not 
transferring  the  z-padding  regions. 

In  our  algorithm  we  use,  instead  of  the  FFTW  transpose,  the  generic  MPI  all-to-all 
communication.  A  potential  drawback  of  this  implementation  is  that  the  MPI  routine  is 
not  executed  in-place,  as  the  FFTW  transpose,  hence  an  additional  storage  is  required. 
However,  the  storage  required  for  the  additional  buffers  is  minor  (see  Eq.(5.7)). 

If  we  use  the  generic  MPI  collective  communication,  an  additional  complication  (but 
leading  to  an  insignificant  cost  increase)  is  the  necessity  of  rearranging  the  FFT  data  before 
transposing  them.  The  entire  procedure  is  illustrated  in  Fig.  13. 

As  in  Fig.  12,  we  start  with  the  FFT  data  distributed  across  the  processors  in  “^-slices” 
(Fig.  13(a)).  The  data  blocks  on  the  individual  processors  have  sizes  (Nx  x  Ny  x  N^), 
i.e. ,  they  are  allocated  with  the  padding  space. 
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Originally,  the  input  data  X  are  stored  only  in  the  “physical”  parts  (Kx  x  Ky  x  K^) 
of  the  blocks,  and  the  remaining  padding  space  is  filled  with  zeros.  After  the  first  FFT,  in 
the  (. x,y )  directions,  is  performed  (by  each  processor  independently),  the  data  expand  to 
(Nx  x  Ny  x  Kz P^),  i.e.,  they  fill  the  padding  in  the  x  and  y  directions,  but  the  2-padding 
remains  unused.  This  is  the  situation  shown  in  Fig.  13(a),  where  the  grey  areas  indicate 
the  actual  data,  and  the  white  areas  the  zero  padding.  The  ordering  of  the  data  elements 
Xn,xiny,nz  (the  nx  changing  the  fastest)  is  shown  by  the  zig-zag  lines  in  the  Figure. 

Fig.  13(b)  shows  a  rearrangement  of  the  data  of  Fig.  13(a),  executed  prior  to  transposi¬ 
tion:  Now,  each  block  of  the  actual  data,  of  size  (Nx  x  Ny  x  K ^),  is  split  into  P  sub-blocks 

of  sizes  (Nx  x  Ny1'1  x  K^),  q  =  0, 1,  ,  . . .  ,  P  —  1.  The  blocks  (marked  “00”  to  “22”)  are 
stored  in  memory  of  the  processor  p  one  after  the  other,  with  the  indicated  element  ordering 
within  each  block  (again,  with  the  fastest  changing  index  nx).  Clearly,  the  rearrangement 
is  carried  out  locally,  and  its  cost  is  proportional  to  Nx  Ny  Kzp>  (no  2-padding  storage  is 
involved). 

Finally,  the  data  layout  after  the  transposition  is  visualized  in  Fig.  13(c).  The  blocks  of 
Fig.  13(b)  have  been  simply  moved  from  processor  to  processor,  with  changing  only  their 
locations  (offsets),  but  not  their  internal  structure  or  ordering.  Due  to  the  offset  changes,  the 
data  blocks  on  each  processors  are  now  “concatenated”  in  the  index  nz,  leaving  contiguous 
zero  padding  in  2  direction.  We  stress  that  only  the  filled  (here,  shaded  in  grey)  blocks, 
and  not  the  empty  padding  space,  are  transferred  between  the  processors. 

The  above  description  applies  to  the  forward  Fourier  transform  of  Eq.(2.12).  In  the 
backward  FFT  (Eq.(2.14))  the  transposition  of  data  is  the  reverse  of  that  in  Fig.  13.  In  this 
case  also  only  the  shaded  blocks  of  Fig.  13(c)  are  transmitted  between  the  processors  to 
obtain  the  configuration  of  Fig.  13(c).  The  data  values  in  the  padding  areas  are  irrelevant, 
and  are  discarded  in  the  process. 
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Figure  13:  Original  FFT  data  (a),  their  rearranged  storage  (b),  and  the  result  of  their 
transposition  by  means  of  the  MPI  all-to-all  communication  routine  (c).  As  before,  the 
zig-zag  lines  indicate  ordering  of  elements. 
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